From statistical modelling to scientific report writing

DSP Handbook for PSB2024

Author

Anna B. Neuheimer

Published

10-11-2024

In this chapter you will learn:
  • how to report your hypothesis testing results (including visualizations) to communicate what your model says about your hypothesis.

1 Reporting your hypothesis testing results

To report your model, you want to communicate a number of things:

  • your best-specified model(s)

  • how well your model explains your response

  • your modelled effects

In a few weeks, we will discuss communicating your hypothesis testing (including results) more formally - i.e. the text and visuals you will use to write your own scientific reports and papers. In the current chapter, we will gather many of the pieces needed for that task.

2 Introducing the examples:

In this chapter, we will report model results for four example best-specified models that showcase a range of different hypotheses.

If you are re-reading this chapter with your own hypothesis you are testing, look for the example that shares the most structure with your own hypothesis.

The examples are:

Example 1: Resp1 ~ Cat1 + 1

Example 2: Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1

Example 3: Resp3 ~ Cont4 + 1

Example 4: Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1

where

  • Resp# is your response variable.
  • Cat# indicates a categorical predictor
  • Cont# indicates a continuous predictor

Note that each example is testing a different hypothesis by fitting a different model to a different data-set. The four examples are not related to one another.

Note also that the examples will all involve a normal error distribution assumption. I will point out when your process would be different for a different error distribution assumption (e.g. Gamma, poisson, and binomial)

Finally, note that I am not going through all the steps that got us to these example best-specified models, but assume we went through the previous steps (Response, Predictors, Hypothesis, Starting model, Model validation and Hypothesis testing) as we practiced in class to choose each best-specified model.

Example 1: Resp1 ~ Cat1 + 1

For example #1, assume your best-specified model says that your response variable (Resp1) is explained by a categorical predictor (Cat1):

Resp1 ~ Cat1

Your best-specified model for example #1 is in an object called bestMod1:

bestMod1

Call:  glm(formula = Resp1 ~ Cat1 + 1, family = gaussian(link = "identity"), 
    data = myDat1)

Coefficients:
(Intercept)      Cat1Sp2      Cat1Sp3  
     20.955        5.877       -4.510  

Degrees of Freedom: 99 Total (i.e. Null);  97 Residual
Null Deviance:      7083 
Residual Deviance: 5210     AIC: 687.1

that was fit to data in myDat1:

str(myDat1)
'data.frame':   100 obs. of  2 variables:
 $ Cat1 : Factor w/ 3 levels "Sp1","Sp2","Sp3": 1 1 1 2 3 3 2 3 2 1 ...
 $ Resp1: num  14.2 11.2 32.9 26.4 20.6 10.9 29 13.6 25.2 29 ...

and the dredge() table you used to pick your bestMod1 is in dredgeOut1:

dredgeOut1
Global model call: glm(formula = Resp1 ~ Cat1 + 1, family = gaussian(link = "identity"), 
    data = myDat1)
---
Model selection table 
  (Intrc) Cat1 df   logLik  AICc delta weight
2   20.95    +  4 -339.547 687.5  0.00      1
1   21.79       2 -354.905 713.9 26.42      0
Models ranked by AICc(x) 

Example 2: Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1

For Example #2, assume your best-specified model says that your response variable (Resp2) is explained by two categorical predictors (Cat2 & Cat3) as well as the interaction between the predictors:

Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1

Your best-specified model for example #2 is in an object called bestMod2:

bestMod2

Call:  glm(formula = Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1, family = gaussian(link = "identity"), 
    data = myDat2)

Coefficients:
          (Intercept)              Cat2TypeB              Cat2TypeC  
               364.91                  37.37                 -75.94  
            Cat2TypeD             Cat3Treat2            Cat3Control  
              -141.97                 -43.63                  64.12  
 Cat2TypeB:Cat3Treat2   Cat2TypeC:Cat3Treat2   Cat2TypeD:Cat3Treat2  
                42.71                  73.61                  66.02  
Cat2TypeB:Cat3Control  Cat2TypeC:Cat3Control  Cat2TypeD:Cat3Control  
                81.25                 -19.08                 -71.72  

Degrees of Freedom: 99 Total (i.e. Null);  88 Residual
Null Deviance:      1249000 
Residual Deviance: 377000   AIC: 1133

that was fit to data in myDat2:

str(myDat2)
'data.frame':   100 obs. of  3 variables:
 $ Resp2: num  278 335 566 341 313 ...
 $ Cat2 : Factor w/ 4 levels "TypeA","TypeB",..: 3 2 2 3 4 1 1 3 4 1 ...
 $ Cat3 : Factor w/ 3 levels "Treat1","Treat2",..: 3 2 3 3 3 1 1 2 1 2 ...

and the dredge() table you used to pick your bestMod2 is in dredgeOut2:

dredgeOut2
Global model call: glm(formula = Resp2 ~ Cat2 + Cat3 + Cat2:Cat3, family = gaussian(link = "identity"), 
    data = myDat2)
---
Model selection table 
  (Int) Ct2 Ct3 Ct2:Ct3     R^2 df   logLik   AICc delta weight
8 364.9   +   +       + 0.69800 13 -553.634 1137.5  0.00  0.964
4 354.9   +   +         0.62540  7 -564.419 1144.1  6.55  0.036
2 388.5   +             0.55300  5 -573.243 1157.1 19.63  0.000
1 348.3                 0.00000  2 -613.508 1231.1 93.64  0.000
3 331.1       +         0.03933  4 -611.502 1231.4 93.92  0.000
Models ranked by AICc(x) 

Example 3: Resp3 ~ Cont4 + 1

For Example #3, assume your best-specified model says that your response variable (Resp3) is explained by one continuous predictor (Cont4):

Resp3 ~ Cont4 + 1

Your best-specified model for example #3 is in an object called bestMod4:

bestMod3

Call:  glm(formula = Resp3 ~ Cont4 + 1, family = gaussian(link = "identity"), 
    data = myDat3)

Coefficients:
(Intercept)        Cont4  
     -226.9        260.1  

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:      276100000 
Residual Deviance: 64930000     AIC: 1628

that was fit to data in myDat4:

str(myDat3)
'data.frame':   100 obs. of  2 variables:
 $ Resp3: num  3668 2774 4594 1094 2100 ...
 $ Cont4: num  13.93 12.15 15.4 4.87 6.64 ...

and the dredge() table you used to pick your bestMod3 is in dredgeOut3:

dredgeOut3
Global model call: glm(formula = Resp3 ~ Cont4 + 1, family = gaussian(link = "identity"), 
    data = myDat3)
---
Model selection table 
  (Intrc) Cont4 df   logLik   AICc  delta weight
2  -226.9 260.1  3 -811.080 1628.4   0.00      1
1  2358.0        2 -883.457 1771.0 142.63      0
Models ranked by AICc(x) 

Example 4: Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1

For Example #4, assume your best-specified model says that your response variable (Resp4) is explained by one categorical predictor (Cat4) and one continuous predictor (Cont6) as well as the interaction between the predictors:

Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1

Your best-specified model for example #3 is in an object called bestMod4:

bestMod4

Call:  glm(formula = Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1, family = gaussian(link = "identity"), 
    data = myDat4)

Coefficients:
    (Intercept)        Cat5Urban         Cat5Wild            Cont6  
      98.346793        -0.236424        -0.867336        -0.002931  
Cat5Urban:Cont6   Cat5Wild:Cont6  
       0.015117         0.030858  

Degrees of Freedom: 99 Total (i.e. Null);  94 Residual
Null Deviance:      4454 
Residual Deviance: 642.8    AIC: 483.9

that was fit to data in myDat4:

str(myDat4)
'data.frame':   100 obs. of  3 variables:
 $ Resp4: num  102 102 104 110 107 ...
 $ Cat5 : Factor w/ 3 levels "Farm","Urban",..: 1 2 2 3 2 1 3 2 3 2 ...
 $ Cont6: num  341 490 643 412 444 ...

and the dredge() table you used to pick your bestMod4 is in dredgeOut4:

dredgeOut4
Global model call: glm(formula = Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1, family = gaussian(link = "identity"), 
    data = myDat4)
---
Model selection table 
   (Int) Ct5       Cn6 Ct5:Cn6    R^2 df   logLik  AICc  delta weight
8  98.35   + -0.002931       + 0.8557  7 -234.930 485.1   0.00      1
4  92.25   +  0.009650         0.8166  5 -246.915 504.5  19.39      0
2  96.93   +                   0.7923  4 -253.133 514.7  29.61      0
3  94.92      0.017780         0.0848  3 -327.278 660.8 175.73      0
1 104.00                       0.0000  2 -331.708 667.5 182.46      0
Models ranked by AICc(x) 

3 Reporting your best-specified models

This section is similar regardless of your model structure (e.g. error distribution assumption). The examples below all assume a normal error distribution assumption but you can use the same process for models of any structure.

Reporting your best-specified model means reporting the terms - the predictors and any interactions - that are in your model. It is good practice to present the model along with the results from the model selection. In this way, you can include multiple best-specified models if there is evidence that more than one might be useful. Depending on your hypothesis and results, you may want to present all models within ∆AICc < 2 of the best-specified model, or all models with any Akaike weight, or simply all models.

Remember from the Hypothesis Testing chapter that you can also use the output from dredge() function to report evidence for how you picked the best-specified model. When we discuss communicating your hypothesis testing in a few weeks, I’ll include some functions that can help make quick work of making a model selection table.

Let’s take a look at how you do this with our examples:

Example 1: Resp1 ~ Cat1 + 1

Global model call: glm(formula = Resp1 ~ Cat1 + 1, family = gaussian(link = "identity"), 
    data = myDat1)
---
Model selection table 
  (Intrc) Cat1 df   logLik  AICc delta weight
2   20.95    +  4 -339.547 687.5  0.00      1
1   21.79       2 -354.905 713.9 26.42      0
Models ranked by AICc(x) 

For Example 1, you will report that your best-specified model is Resp1 ~ Cat1 + 1, i.e. that there is evidence that Cat1 explains variability in Resp1. This was chosen as the best-specified model as it had the lowest AICc and the next highest rank model had a ∆AICc of 26.42 (i.e. ∆AICc > 2). The best-specified model had an Akaike weight of 1.

Example 2: Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1

Global model call: glm(formula = Resp2 ~ Cat2 + Cat3 + Cat2:Cat3, family = gaussian(link = "identity"), 
    data = myDat2)
---
Model selection table 
  (Int) Ct2 Ct3 Ct2:Ct3     R^2 df   logLik   AICc delta weight
8 364.9   +   +       + 0.69800 13 -553.634 1137.5  0.00  0.964
4 354.9   +   +         0.62540  7 -564.419 1144.1  6.55  0.036
2 388.5   +             0.55300  5 -573.243 1157.1 19.63  0.000
1 348.3                 0.00000  2 -613.508 1231.1 93.64  0.000
3 331.1       +         0.03933  4 -611.502 1231.4 93.92  0.000
Models ranked by AICc(x) 

For Example 2, you will report that your best-specified model is Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1, i.e. that there is evidence that Cat2 and Cat3 explain variability in Resp2, and that there is an interaction between the effect of Cat2 and Cat3 on your response - i.e. the effect of Cat2 on Resp2 depends on the value of Cat3.

This was chosen as the best-specified model as it had the lowest AICc and the next highest rank model had a ∆AICc of 6.55 (i.e. ∆AICc > 2). The best-specified model had an Akaike weight of 0.964.

Example 3: Resp3 ~ Cont4 + 1

Global model call: glm(formula = Resp3 ~ Cont4 + 1, family = gaussian(link = "identity"), 
    data = myDat3)
---
Model selection table 
  (Intrc) Cont4 df   logLik   AICc  delta weight
2  -226.9 260.1  3 -811.080 1628.4   0.00      1
1  2358.0        2 -883.457 1771.0 142.63      0
Models ranked by AICc(x) 

For Example 3, you will report that your best-specified model is Resp3 ~ Cont4 + 1, i.e. that there is evidence that Cont4 explains variability in Resp3.

This was chosen as the best-specified model as it had the lowest AICc and the next highest rank model had a ∆AICc of 142.63 (i.e. ∆AICc > 2). The best-specified model had an Akaike weight of 1.

Example 4: Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1

Global model call: glm(formula = Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1, family = gaussian(link = "identity"), 
    data = myDat4)
---
Model selection table 
   (Int) Ct5       Cn6 Ct5:Cn6    R^2 df   logLik  AICc  delta weight
8  98.35   + -0.002931       + 0.8557  7 -234.930 485.1   0.00      1
4  92.25   +  0.009650         0.8166  5 -246.915 504.5  19.39      0
2  96.93   +                   0.7923  4 -253.133 514.7  29.61      0
3  94.92      0.017780         0.0848  3 -327.278 660.8 175.73      0
1 104.00                       0.0000  2 -331.708 667.5 182.46      0
Models ranked by AICc(x) 

For Example 4, you will report that your best-specified model is Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1. This model says that there is evidence that Cat5 and Cont6 explain variability in Resp4, and that there is an interaction between Cat5 and Cont6 - i.e. the effect of Cont6 on Resp4 depends on the value of Cat5.

This was chosen as the best-specified model as it had the lowest AICc and the next highest rank model had a ∆AICc of 19.39 (i.e. ∆AICc > 2). The best-specified model had an Akaike weight of 1.

4 Reporting how well your model explains your response

This section is similar regardless of your model structure (e.g. error distribution assumption). The examples below all assume a normal error distribution assumption but you can use the same process for models of any structure.

In this section you will determine

  • how much deviance in your response is explained by your model

  • how important each of your predictors is in explaining that deviance

4.1 how much variation in your response is explained by your model

If you will recall, your whole motivation for pursuing statistical modelling was to explain variation in your response. Thus, it is important that you quantify how much variation in your response you are able to explain by your model.

Note that we will discuss this as “explained deviance” rather than “explained variation”. This is because the term “variance” is associated with models where the error distribution assumption is normal, whereas deviance is a more universal term.

When you have a normal error distribution assumption and linear shape assumption, you can capture the amount of explained deviance simply by comparing the variance in your response (i.e. the starting variation before you fit your model - the null variation) vs. the variance in your model residuals (i.e. the remaining variation after you fit your model - the residual variation) as the \(R^2\):

\(R^2 = 1 - \frac{residual variation}{null variation}\)

From this equation, you can see how, if your model is able to explain all the variation in the response, the residual variation will be zero, and \(R^2 = 1\). Alternatively, if the model explains no variation in the response the residual variation equals the null variation and \(R^2 = 0\).

For models with other error distribution and shape assumptions, you need another way of estimating the goodness of fit of your model. You can do this through a pseudo \(R^2\).

One useful pseudo \(R^2\) is called the Likelihood Ratio \(R^2\) or \(R^2_{LR}\). The \(R^2_{LR}\) compares the likelihood of your best-specified model to the likelihood of the null model:

\(R^2_{LR} = 1-exp(-\frac{2}{n}(log𝓛(model)-log𝓛(null)))\)

where \(n\) is the number of observations, \(log𝓛(model)\) is the log-likelihood of your model, and \(log𝓛(null)\) is the log-likelihood of the null model. The \(R^2_{LR}\) is the type of pseudo \(R^2\) that shows up in your dredge() output when you add extra = "R^2" to the dredge() call. You can calculate \(R^2_{LR}\) by hand, read it from our dredge() output, or estimate it using r.squaredLR() from the MuMIn package:

r.squaredLR(bestMod1)
[1] 0.2644618
attr(,"adj.r.squared")
[1] 0.2646806

Note here that two values of \(R^2_{LR}\) are reported. The adjusted pseudo \(R^2\) given here under attr(,"adj.r.squared") has been scaled so that \(R^2_{LR}\) can reach a maximum of 1, similar to a regular \(R^2\)1.

Let’s compare this to the traditional \(R^2\):

1-summary(bestMod1)$deviance/summary(bestMod1)$null.deviance
[1] 0.2644618

Note that the two estimates are similar: One nice feature of the \(R^2_{LR}\) is that it is equivalent to the regular \(R^2\) when our model assumes a normal error distribution assumption and linear shape assumption, so you can use \(R^2_{LR}\) for any of the models that we are discussing in class.

4.2 how important each of your predictors is in explaining that variation

When you have more than one predictor in your model. you may also want to report how relatively important each predictor is to explaining deviance in your response. This is also called “partitioning the explained deviance” to the predictors.

To get an estimate of how much deviance in your response one particular predictor explains, you may be tempted to compute the explained deviance (\(R^2\)) estimates of models fit to data with and without that particular predictor. Let’s try with our Example 4:

dredgeOut4
Global model call: glm(formula = Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1, family = gaussian(link = "identity"), 
    data = myDat4)
---
Model selection table 
   (Int) Ct5       Cn6 Ct5:Cn6    R^2 df   logLik  AICc  delta weight
8  98.35   + -0.002931       + 0.8557  7 -234.930 485.1   0.00      1
4  92.25   +  0.009650         0.8166  5 -246.915 504.5  19.39      0
2  96.93   +                   0.7923  4 -253.133 514.7  29.61      0
3  94.92      0.017780         0.0848  3 -327.278 660.8 175.73      0
1 104.00                       0.0000  2 -331.708 667.5 182.46      0
Models ranked by AICc(x) 

If you want to get an estimate as to how much response deviance the Cont6 predictor explains, you might compare the \(R^2\) of a model with and without the Cont6 predictor.

Let’s compare comparing model #4 (that includes Cont6) and model #2 (that doesn’t include Cont6):

R2.mod4 <- (dredgeOut4$`R^2`[2]) # model #4 is the second row in dredgeOut4
R2.mod2 <- (dredgeOut4$`R^2`[3]) # model #2 is the third row in dredgeOut4

diffR2 <- R2.mod4 - R2.mod2 # find estimated contribution of Cont to explained deviance

diffR2
[1] 0.02429225

So it looks like 2.4% of the variability in Resp4 is explained by Cont6.

But what if instead you chose to compare model #3 (that includes Cont6) and model #1 (that doesn’t include Cont6):

R2.mod3 <- (dredgeOut4$`R^2`[4]) # model #3 is the fourth row in dredgeOut4
R2.mod1 <- (dredgeOut4$`R^2`[5]) # model #1 is the fifth row in dredgeOut4

diffR2 <- R2.mod3 - R2.mod1 # find estimated contribution of Cont to explained deviance

diffR2
[1] 0.08480009

Now it looks like 8.5% of the variability in Resp4 is explained by Cont6! Quite a different answer! Your estimates of the contribution of Cont6 to explaining the response deviation don’t agree because of collinearity among our predictors2.

There are a few options you can use to get a robust estimate of how much each predictor is contributing to explained deviance in your response.

One option for partitioning the explained deviance when you have collinearity among your predictors is hierarchical partitioning. Hierarchical partitioning estimates the average independent contribution of each predictor to the total explained variance by considering all models in the candidate model set3. This method is beyond the scope of the course but see the rdacca.hp package for an example of how to do this.

Another method (that we will be using) for estimating the importance of each term (predictor or interaction) in your model is by again looking at the candidate model set ranking made by dredge(). Here you can measure the importance of a predictor by summing up the Akaike weights for any model that includes a particular predictor. The Akaike weight of a model compares the likelihood of the model scaled to the total likelihood of all the models in the candidate model set. The sum of Akaike weights for models including a particular predictor tells you how important the predictor is in explaining the deviance in your response. You can calculate the sum of Akaike weights directly with the sw() function in the MuMIn package:

sw(dredgeOut4)
                     Cat5 Cont6 Cat5:Cont6
Sum of weights:      1    1     1         
N containing models: 3    3     1         

Here we can see that all model terms (the predictors Cat5 and Cont6 as well as the interaction) are equally important in explaining the deviance in Resp4 (they appear in all models that have any Akaike weight).

Let’s look at these two steps

  • how much deviance in your response is explained by your model

  • how important each of your predictors is in explaining that deviance

with our examples:

Example 1: Resp1 ~ Cat1 + 1

  • how much deviance in your response is explained by your model
R2 <- r.squaredLR(bestMod1)

R2
[1] 0.2644618
attr(,"adj.r.squared")
[1] 0.2646806

The best-specified model explains 26.4% of the deviance in Resp1.

  • how important each of your predictors is in explaining that deviance
dredgeOut2
Global model call: glm(formula = Resp2 ~ Cat2 + Cat3 + Cat2:Cat3, family = gaussian(link = "identity"), 
    data = myDat2)
---
Model selection table 
  (Int) Ct2 Ct3 Ct2:Ct3     R^2 df   logLik   AICc delta weight
8 364.9   +   +       + 0.69800 13 -553.634 1137.5  0.00  0.964
4 354.9   +   +         0.62540  7 -564.419 1144.1  6.55  0.036
2 388.5   +             0.55300  5 -573.243 1157.1 19.63  0.000
1 348.3                 0.00000  2 -613.508 1231.1 93.64  0.000
3 331.1       +         0.03933  4 -611.502 1231.4 93.92  0.000
Models ranked by AICc(x) 
sw(dredgeOut1)
                     Cat1
Sum of weights:      1   
N containing models: 1   

With Example 1, you have only one predictor (Cat1) and so this predictor is responsible for explaining all of the variability in your response (Resp1). You can see that it appears in all models with any weight with your sw() function from the MuMIn package.

Example 2: Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1

  • how much deviance in your response is explained by your model
R2 <- r.squaredLR(bestMod2)

R2
[1] 0.6980474
attr(,"adj.r.squared")
[1] 0.6980507

The best-specified model explains 69.8% of the deviance in Resp2.

  • how important each of your predictors is in explaining that deviance
dredgeOut2
Global model call: glm(formula = Resp2 ~ Cat2 + Cat3 + Cat2:Cat3, family = gaussian(link = "identity"), 
    data = myDat2)
---
Model selection table 
  (Int) Ct2 Ct3 Ct2:Ct3     R^2 df   logLik   AICc delta weight
8 364.9   +   +       + 0.69800 13 -553.634 1137.5  0.00  0.964
4 354.9   +   +         0.62540  7 -564.419 1144.1  6.55  0.036
2 388.5   +             0.55300  5 -573.243 1157.1 19.63  0.000
1 348.3                 0.00000  2 -613.508 1231.1 93.64  0.000
3 331.1       +         0.03933  4 -611.502 1231.4 93.92  0.000
Models ranked by AICc(x) 
sw(dredgeOut2)
                     Cat2 Cat3 Cat2:Cat3
Sum of weights:      1.00 1.00 0.96     
N containing models:    3    3    1     

Here you can see that Cat2 and Cat3 are equally important in explaining the deviance in Resp2 (they appear in all models that have any Akaike weight), while the interaction term between Cat2 and Cat3 is less important (only appearing in one model with Akaike weight, though this is the top model).

Example 3: Resp3 ~ Cont4 + 1

  • how much deviance in your response is explained by your model
R2 <- r.squaredLR(bestMod3)

R2
[1] 0.764852
attr(,"adj.r.squared")
[1] 0.764852

The best-specified model explains 76.5% of the deviance in Resp3.

  • how important each of your predictors is in explaining that deviance
dredgeOut3
Global model call: glm(formula = Resp3 ~ Cont4 + 1, family = gaussian(link = "identity"), 
    data = myDat3)
---
Model selection table 
  (Intrc) Cont4 df   logLik   AICc  delta weight
2  -226.9 260.1  3 -811.080 1628.4   0.00      1
1  2358.0        2 -883.457 1771.0 142.63      0
Models ranked by AICc(x) 
sw(dredgeOut3)
                     Cont4
Sum of weights:      1    
N containing models: 1    

With Example 3, you have only one predictor (Cont6) and so this predictor is responsible for explaining all of the variability in your response (Resp3). You can see that it appears in all models with any weight with your sw() function from the MuMIn package.

Example 4: Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1

  • how much deviance in your response is explained by your model
R2 <- r.squaredLR(bestMod4)

R2
[1] 0.855657
attr(,"adj.r.squared")
[1] 0.8567834

The best-specified model explains 85.6% of the deviance in Resp2.

  • how important each of your predictors is in explaining that deviance
dredgeOut4
Global model call: glm(formula = Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1, family = gaussian(link = "identity"), 
    data = myDat4)
---
Model selection table 
   (Int) Ct5       Cn6 Ct5:Cn6    R^2 df   logLik  AICc  delta weight
8  98.35   + -0.002931       + 0.8557  7 -234.930 485.1   0.00      1
4  92.25   +  0.009650         0.8166  5 -246.915 504.5  19.39      0
2  96.93   +                   0.7923  4 -253.133 514.7  29.61      0
3  94.92      0.017780         0.0848  3 -327.278 660.8 175.73      0
1 104.00                       0.0000  2 -331.708 667.5 182.46      0
Models ranked by AICc(x) 
sw(dredgeOut4)
                     Cat5 Cont6 Cat5:Cont6
Sum of weights:      1    1     1         
N containing models: 3    3     1         

Here you can see that Cat5, Cont6 and the interaction Cat5:Cont6 are all equally important in explaining the deviance in Resp4 (they appear in all models that have any Akaike weight).

5 Reporting your modelled effects

In the last chapter of these notes, we discussed testing your hypothesis.

Remember that testing your hypothesis is looking for evidence of effects of your predictors on your response. By testing your hypothesis, you determined that there was evidence that the effects of your predictors on your response are non-zero.

Here you will gather the information to communicate visually and quantitatively the magnitude and direction of your effects and determine (in some cases) where effects are similar or different from one another, answering questions like:

  • how much does a change (effect) in your predictor change your response?

  • is that change (effect) positive (your response) or negative (your response decreases)?

  • is that change (effect) similar across levels of your categorical predictor?

You will report your modelled effects in two ways:

  1. Visualizing your model effects: Here, you will make plots of your model effects including the effects themselves, uncertainty around the effects and your observations. These plots will help communicate your modelled effects as well as how well your model fits your observations.

  2. Quantifying your model effects: Here, you will report (in numbers) the magnitude and direction of your modelled effects along with the uncertainty. Exactly how you do that will depend on the structure of your model (e.g. the error distribution assumption). We will go through a number of examples here. You will also look for evidence of whether modelled effects of a categorical predictor differ across all levels.

5.1 Visualizing your model effects

This section is similar regardless of your model structure (e.g. error distribution assumption). The examples below all assume a normal error distribution assumption but you can use the same process for models of any structure.

Here, you will visualize how each predictor is affecting the response by drawing how your modelled fitted values (the estimates of your response made by the model) change with your predictors, along with a measure of uncertainty around those fitted values. You also want to include your data (actual observations of your response and predictors) so you can start to communicate how well your model captures your hypothesis, and what amount of deviance in your response is explained by your model.

There are a lot of different R packages that make it easy to quickly visualize your model. We will focus on two methods that will allow you to make quick visualizations of your model and/or customize the figures as you would like.

  • visualizing using the visreg package: This will let you get a quick look at your model object with a limited amount of customization.

  • vizualizing “by hand”: Here, “by hand” is a bit of a silly description as R will be doing the work for you. What I mean by “by hand” is that you will be building the plot yourself by querying your model object. This method is very flexible and will lead you to a fully customizable visualization. This process involves i) choosing the values of your predictors at which to make estimates of your fitted values, ii) using predict() to use your model to estimate your response variable at those values of your predictors, and iii) use the model estimates to plot your model fit. Aside: this is similar to what we will be doing when we learn about model predictions in the next chapter of these notes.

Below are visualizations for each of our examples.

Example 1: Resp1 ~ Cat1 + 1

Visualizing with the visreg package

library (visreg) # load visreg package for visualization

library(ggplot2) # load ggplot2 package for visualization

visreg(bestMod1, # model to visualize
       scale = "response", # plot on the scale of the response
       xvar = "Cat1", # predictor on x-axis
       #by = ..., # if you want to include a 2nd predictor plotted as colour
       #breaks = ..., # if you want to control how the colour predictor is plotted
       #cond = , # if you want to include a 3rd predictor
       #overlay = TRUE, # to plot as overlay or panels, when there is 
       #rug = FALSE, # to turn off the rug. The rug shows you where you have positive (appearing on the top axis) and negative (appearing on the bottom axis) residuals
       gg = TRUE)+ # to plot as a ggplot, with ggplot, you will have more control over the look of the plot.
  ylab("Response, (units)")+ # change y-axis label
  xlab("Cat1, (units)")+ # change x-axis label
  theme_bw() # change ggplot theme

Unfortunately due to limitations of the visreg package, you can not easily add you observations onto plots where the x-axis is a categorical predictor. But that’s ok, because there are other options…

Visualizing by hand

To plot by hand, you will

  1. first make a data frame containing the value of your predictors at which you want to plot effects on the response:
# Set up your predictors for the visualized fit
forCat1<-unique(myDat1$Cat1) # every value of your categorical predictor

# create a data frame with your predictors
forVis<-expand.grid(Cat1=forCat1) # expand.grid() function makes sure you have all combinations of predictors
  1. Next, you will use the predict() function4 to get the model estimates of your response variable at those values of your predictors:
# Get your model fit estimates at each value of your predictors
modFit<-predict(bestMod1, # the model
                newdata = forVis, # the predictor values
                type = "response", # make the predictions on the response scale
                se.fit = TRUE) # include uncertainty estimate

forVis$Fit<-modFit$fit # add your fit to the data frame
forVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frame
forVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame
  1. Finally, you will use these model estimates to make your plot:
library(ggplot2)

ggplot() + # start ggplot
  geom_point(data = myDat1, # add observations to your plot
             mapping = aes(x = Cat1, y = Resp1), 
             position=position_jitter(width=0.1)) + # control position of data points so they are easier to see on the plot
  geom_errorbar(data = forVis, # add the uncertainty to your plot
              mapping = aes(x = Cat1, y = Fit, ymin = Lower, ymax = Upper),
              linewidth=1.2) + # control thickness of errorbar line
  geom_point(data = forVis, # add the modelled fit to your plot
              mapping = aes(x = Cat1, y = Fit), 
             shape = 22, # shape of point
             size = 3, # size of point
             fill = "white", # fill colour for plot
             col = 'black') + # outline colour for plot
  ylab("Resp1, (units)")+ # change y-axis label
  xlab("Cat1, (units)")+ # change x-axis label
  theme_bw() # change ggplot theme

Example 2: Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1

Visualizing with the visreg package

As you have more than one predictor in Example 2, there are a number of ways you can visualize your effects.

Here is an example of a visualization with Cat2 on the x-axis and the effects of Cat3 plotted as separate panels:

library (visreg) # load visreg package for visualization

library(ggplot2) # load ggplot2 package for visualization

visreg(bestMod2, # model to visualize
       scale = "response", # plot on the scale of the response
       xvar = "Cat2", # predictor on x-axis
       by = "Cat3", # if you want to include a 2nd predictor plotted as colour
       #breaks = ..., # if you want to control how the colour predictor is plotted
       #cond = , # if you want to include a 3rd predictor
       overlay = FALSE, # to plot as overlay or panels, when there is 
       #rug = FALSE, # to turn off the rug. The rug shows you where you have positive (appearing on the top axis) and negative (appearing on the bottom axis) residuals
       gg = TRUE)+ # to plot as a ggplot, with ggplot, you will have more control over the look of the plot.
  ylab("Response, (units)")+ # change y-axis label
  xlab("Cat2, (units)")+ # change x-axis label
  theme_bw() # change ggplot theme

# tip: we can't include our observations on this plot due to the limitations of the visreg package when plotting with a categorical predictor on the x-axis

And here is an example of a visualization with Cat2 on the x-axis and the effects of Cat3 overlayed on the plot as different colours for each category (level) of Cat3:

visreg(bestMod2, # model to visualize
       scale = "response", # plot on the scale of the response
       xvar = "Cat2", # predictor on x-axis
       by = "Cat3", # if you want to include a 2nd predictor plotted as colour
       #breaks = ..., # if you want to control how the colour predictor is plotted
       #cond = , # if you want to include a 3rd predictor
       overlay = TRUE, # to plot as overlay or panels, when there is 
       #rug = FALSE, # to turn off the rug. The rug shows you where you have positive (appearing on the top axis) and negative (appearing on the bottom axis) residuals
       gg = TRUE)+ # to plot as a ggplot, with ggplot, you will have more control over the look of the plot.
  ylab("Response, (units)")+ # change y-axis label
  xlab("Cat2, (units)")+ # change x-axis label
  theme_bw() # change ggplot theme

# tip: we can't include our observations on this plot due to the limitations of the visreg package when plotting with a categorical predictor on the x-axis

Unfortunately due to limitations of the visreg package, you can not easily add you observations onto plots where the x-axis is a categorical predictor. But that’s ok, because there are other options…

Visualizing by hand

#### i) choosing the values of your predictors at which to make a prediction

# Set up your predictors for the visualized fit
forCat2<-unique(myDat2$Cat2) # every level of your categorical predictor
forCat3<-unique(myDat2$Cat3) # every level of your categorical predictor
  
# create a data frame with your predictors
forVis<-expand.grid(Cat2 = forCat2, Cat3 = forCat3) # expand.grid() function makes sure you have all combinations of predictors

#### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor


# Get your model fit estimates at each value of your predictors
modFit<-predict(bestMod2, # the model
                newdata = forVis, # the predictor values
                type = "response", # make the predictions on the response scale
                se.fit = TRUE) # include uncertainty estimate

forVis$Fit<-modFit$fit # add your fit to the data frame
forVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frame
forVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame

#### iii) use the model estimates to plot your model fit

library(ggplot2) # load ggplot2 library

ggplot() + # start ggplot
  geom_point(data = myDat2, # add observations to your plot
             mapping = aes(x = Cat2, y = Resp2, col = Cat3), 
             position=position_jitterdodge(jitter.width=0.75, dodge.width=0.75)) + # control position of data points so they are easier to see on the plot
  geom_errorbar(data = forVis, # add the uncertainty to your plot
              mapping = aes(x = Cat2, y = Fit, ymin = Lower, ymax = Upper, col = Cat3),
              position=position_dodge(width=0.75), # control position of data points so they are easier to see on the plot
              size=1.2) + # control thickness of errorbar line
  geom_point(data = forVis, # add the modelled fit to your plot
              mapping = aes(x = Cat2, y = Fit, fill = Cat3), 
             shape = 22, # shape of point
             size=3, # size of point
             col = 'black', # outline colour for point
             position=position_dodge(width=0.75)) + # control position of data points so they are easier to see on the plot
  ylab("Resp2, (units)")+ # change y-axis label
  xlab("Cat2, (units)")+ # change x-axis label
  labs(fill="Cat3, (units)", col="Cat3, (units)") + # change legend title
  theme_bw() # change ggplot theme

Example 3: Resp3 ~ Cont4 + 1

Visualizing with the visreg package

Notice how you can include the gg = TRUE argument to plot this as a ggplot type figure. This allows you to add your data onto the visualization of your model.

library (visreg) # load visreg package for visualization

library(ggplot2) # load ggplot2 package for visualization

visreg(bestMod3, # model to visualize
       scale = "response", # plot on the scale of the response
       xvar = "Cont4", # predictor on x-axis
       #by = ..., # if you want to include a 2nd predictor plotted as colour
       #breaks = ..., # if you want to control how the colour predictor is plotted
       #cond = , # if you want to include a 3rd predictor
       #overlay = TRUE, # to plot as overlay or panels, when there is 
       #rug = FALSE, # to turn off the rug. The rug shows you where you have positive (appearing on the top axis) and negative (appearing on the bottom axis) residuals
       gg = TRUE)+ # to plot as a ggplot, with ggplot, you will have more control over the look of the plot.
  geom_point(data = myDat3,
             mapping = aes(x = Cont4, y = Resp3))+
  ylab("Response, (units)")+ # change y-axis label
  xlab("Cont4, (units)")+ # change x-axis label
  theme_bw() # change ggplot theme

# Note: you CAN include your observations on your plot with visreg when you have a continuous predictor on the x-axis

Visualizing by hand

#### i) choosing the values of your predictors at which to make a prediction


# Set up your predictors for the visualized fit
forCont4<-seq(from = min(myDat3$Cont4), to = max(myDat3$Cont4), by = 1)# a sequence of numbers in your continuous predictor range
  
# create a data frame with your predictors
forVis<-expand.grid(Cont4 = forCont4) # expand.grid() function makes sure you have all combinations of predictors.  

#### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor


# Get your model fit estimates at each value of your predictors
modFit<-predict(bestMod3, # the model
                newdata = forVis, # the predictor values
                type = "response", # make the predictions on the response scale
                se.fit = TRUE) # include uncertainty estimate

forVis$Fit<-modFit$fit # add your fit to the data frame
forVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frame
forVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame

#### iii) use the model estimates to plot your model fit


library(ggplot2) # load ggplot2 library

ggplot() + # start ggplot
  
  geom_point(data = myDat3, # add observations to your plot
             mapping = aes(x = Cont4, y = Resp3)) + # control position of data points so they are easier to see on the plot
  
  geom_line(data = forVis, # add the modelled fit to your plot
              mapping = aes(x = Cont4, y = Fit),
              size = 1.2) + # control thickness of line
  
    geom_line(data = forVis, # add uncertainty to your plot (upper line)
              mapping = aes(x = Cont4, y = Upper),
              size = 0.8, # control thickness of line
              linetype = 2) + # control style of line
  
      geom_line(data = forVis, # add uncertainty to your plot (lower line)
              mapping = aes(x = Cont4, y = Lower),
              size = 0.8, # control thickness of line
              linetype = 2) + # control style of line
  
  ylab("Resp3, (units)") + # change y-axis label
  
  xlab("Cont4, (units)") + # change x-axis label
  
  theme_bw() # change ggplot theme

Example 4: Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1

Visualizing with the visreg package

As you will see in @sec_quantifying, there is no evidence of an effect of Cont6 on Resp3 when Cat5 = Farm - i.e. the coefficient (slope) is not different from zero. In such cases, you would typically not include the effect on the plot (red line in the next figure), but I include it here to illustrate the visualization techniques.

A plot with Cont6 on the x-axis:

library (visreg) # load visreg package for visualization

library(ggplot2) # load ggplot2 package for visualization

visreg(bestMod, # model to visualize
       scale = "response", # plot on the scale of the response
       xvar = "Cont6", # predictor on x-axis
       by = "Cat5", # predictor plotted as colour
       #breaks = 3, # if you want to control how the colour predictor is plotted
       #cond = , # if you want to include a 4th predictor
       overlay = TRUE, # to plot as overlay or panels 
       rug = FALSE, # to include a rug
       gg = TRUE)+ # to plot as a ggplot
  geom_point(data = myDat4, # data
             mapping = aes(x = Cont6, y = Resp4, col = Cat5))+ # add data to your plot
  #ylim(0, 60)+ # adjust the y-axis units
  ylab("Response, (units)")+ # change y-axis label
  xlab("Cont6, (units)")+ # change x-axis label
  theme_bw() # change ggplot theme

A plot with Cat5 on the x-axis: Notice how the continuous predictor is represented by different levels in the colours on the plot. Here you’ve asked for three levels with breaks = 3. Note that you can not include your observations on the visreg plot when the x-axis predictor is a category

visreg(bestMod, # model to visualize
       scale = "response", # plot on the scale of the response
       xvar = "Cat5", # predictor on x-axis
       by = "Cont6", # predictor plotted as colour
       breaks = 3, # if you want to control how the colour predictor is plotted
       #cond = , # if you want to include a 4th predictor
       overlay = TRUE, # to plot as overlay or panels 
       rug = FALSE, # to include a rug
       gg = TRUE)+ # to plot as a ggplot
  ylab("Response, (units)")+ # change y-axis label
  xlab("Cat5, (units)")+ # change x-axis label
  labs(fill="Cont6, (units)", col="Cont6, (units)") + # change legend title
  theme_bw() # change ggplot theme

You can also specify at which levels the breaks should occur with the breaks = ... argument. For example, you can ask visreg to plot the modelled effects when Cont6 = 400 and Cont6 = 600:

visreg(bestMod, # model to visualize
       scale = "response", # plot on the scale of the response
       xvar = "Cat5", # predictor on x-axis
       by = "Cont6", # predictor plotted as colour
       breaks = c(400,600), # if you want to control how the colour predictor is plotted
       #cond = , # if you want to include a 4th predictor
       overlay = TRUE, # to plot as overlay or panels 
       rug = FALSE, # to include a rug
       gg = TRUE)+ # to plot as a ggplot
  ylab("Response, (units)")+ # change y-axis label
  xlab("Cat5, (units)")+ # change x-axis label
  labs(fill="Cont6, (units)", col="Cont6, (units)") + # change legend title
  theme_bw() # change ggplot theme

Note that you can not include your observations on the visreg plot when the x-axis predictor is a category.

Visualizing by hand

Here’s an example with Cont6 on the x-axis:

#### i) choosing the values of your predictors at which to make a prediction


# Set up your predictors for the visualized fit
forCat5<-unique(myDat4$Cat5) # all levels of your categorical predictor
forCont6<-seq(from = min(myDat4$Cont6), to = max(myDat4$Cont6), by = 1)# a sequence of numbers in your continuous predictor range
  
# create a data frame with your predictors
forVis<-expand.grid(Cat5 = forCat5, Cont6 = forCont6) # expand.grid() function makes sure you have all combinations of predictors.  

#### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor


# Get your model fit estimates at each value of your predictors
modFit<-predict(bestMod4, # the model
                newdata = forVis, # the predictor values
                type = "response", # make the predictions on the response scale
                se.fit = TRUE) # include uncertainty estimate

forVis$Fit<-modFit$fit # add your fit to the data frame
forVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frame
forVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame

#### iii) use the model estimates to plot your model fit


library(ggplot2) # load ggplot2 library

ggplot() + # start ggplot
  
  geom_point(data = myDat4, # add observations to your plot
             mapping = aes(x = Cont6, y = Resp4, col = Cat5)) + # control position of data points so they are easier to see on the plot
  
  geom_line(data = forVis, # add the modelled fit to your plot
              mapping = aes(x = Cont6, y = Fit, col = Cat5),
              size = 1.2) + # control thickness of line
  
    geom_line(data = forVis, # add uncertainty to your plot (upper line)
              mapping = aes(x = Cont6, y = Upper, col = Cat5),
              size = 0.4, # control thickness of line
              linetype = 2) + # control style of line
  
      geom_line(data = forVis, # add uncertainty to your plot (lower line)
              mapping = aes(x = Cont6, y = Lower, col = Cat5),
              size = 0.4, # control thickness of line
              linetype = 2) + # control style of line
  
  ylab("Resp4, (units)") + # change y-axis label
  
  xlab("Cont6, (units)") + # change x-axis label
  
  labs(fill="Cat5, (units)", col="Cat5, (units)") + # change legend title
  
  theme_bw() # change ggplot theme

Here’s an example with Cat5 on the x-axis:

#### i) choosing the values of your predictors at which to make a prediction

# Set up your predictors for the visualized fit
forCat5<-unique(myDat4$Cat5) # all levels of your categorical predictor
forCont6<-c(400, 600) # a sequence of numbers in your continuous predictor range
  
# create a data frame with your predictors
forVis<-expand.grid(Cat5 = forCat5, Cont6 = forCont6) # expand.grid() function makes sure you have all combinations of predictors.  

#### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor

# Get your model fit estimates at each value of your predictors
modFit<-predict(bestMod4, # the model
                newdata = forVis, # the predictor values
                type = "response", # make the predictions on the response scale
                se.fit = TRUE) # include uncertainty estimate

forVis$Fit<-modFit$fit # add your fit to the data frame
forVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frame
forVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame

#### iii) use the model estimates to plot your model fit


library(ggplot2) # load ggplot2 library


ggplot() + # start ggplot
  
  geom_point(data = myDat4, # add observations to your plot
             mapping = aes(x = Cat5, y = Resp4, col = Cont6), 
             position=position_jitterdodge(jitter.width=0.75, dodge.width=0.75)) + # control position of data points so they are easier to see on the plot
  
  geom_errorbar(data = forVis, # add the uncertainty to your plot
              mapping = aes(x = Cat5, y = Fit, ymin = Lower, ymax = Upper, col = Cont6),
              position=position_dodge(width=0.75), # control position of data points so they are easier to see on the plot
              size=1.2) + # control thickness of errorbar line
  
  geom_point(data = forVis, # add the modelled fit to your plot
              mapping = aes(x = Cat5, y = Fit, fill = Cont6), 
             shape = 22, # shape of point
             size=3, # size of point
             col = 'black', # outline colour for point
             position=position_dodge(width=0.75)) + # control position of data points so they are easier to see on the plot
  
  ylab("Resp4, (units)")+ # change y-axis label
  
  xlab("Cat5, (units)")+ # change x-axis label
  
  labs(fill="Cont6, (units)", col="Cont6, (units)") + # change legend title
  
  theme_bw() # change ggplot theme

You might notice that many of the plots above are communicating 3-dimensions (one response + two predictors) in a 2-dimensional plot. There are other ways of making 3-dimensional plots in R, e.g. with the visreg package using the visreg2d() function in the visreg package:

visreg2d(fit = bestMod4, # your model
         xvar = "Cont6", # what to plot on the x-axis
         yvar = "Cat5", # what to plot on the y-axis
         scale = "response") # make sure fitted values (colours) are on the scale of the response

or “by hand” using the geom_raster() function in the ggplot2 package:

# Set up your predictors for the visualized fit
forCont6<-seq(from = min(myDat4$Cont6), 
             to = max(myDat4$Cont6), 
             by = 0.1) # e.g. a sequence of Cont values
forCat5<-unique(myDat4$Cat5) # every value of your categorical predictor

# create a data frame with your predictors
forVis<-expand.grid(Cont6=forCont6, Cat5=forCat5) # expand.grid() function makes sure you have all combinations of predictors


# Get your model fit estimates at each value of your predictors
modFit<-predict(bestMod4, # the model
                newdata = forVis, # the predictor values
                type = "response", # make the predictions on the response scale
                se.fit = TRUE) # include uncertainty estimate

forVis$Fit<-modFit$fit # add your fit to the data frame
forVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frame
forVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame

# create your plot
ggplot() + # start your ggplot
  geom_raster(data = forVis, aes(x = Cont6, y = Cat5, fill = Fit))+ # add the 3 dimensions as a raster
  geom_point(data = myDat4, aes(x = Cont6, y = Cat5, colour = Resp)) # add your data

EXTRA: As you can see, these plots represent the 3rd dimension by using colour. You can also make actual 3 dimensional plots in R with the plotly package. These plots are interactive which makes them more useful than static 3d plots. Click on the plot and move your mouse to rotate the plot!

library(plotly) # load the plotly package

plot_ly(data = forVis, # the data with your model predictions (made above)
        x = ~Cont6, # what to plot on the x axis
        y = ~Cat5, # what to plot on the y axis
        z = ~Fit, # what to plot on the z axis
        type="scatter3d", # type of plot
        mode="markers") %>% # type of plot
  add_markers(data = myDat4, x = ~Cont6, y = ~Cat5, z = ~Resp4) # add your data

5.2 Quantifying your model effects

How you quantify your model effects varies with your model structure. If it is your first time reading this, read through the examples that present models assuming a normal error distribution assumption. This will help you understand why we are reporting modelled effects in this way. Then, (if relevant) look at the section on the error distribution assumption you are interested in for your model.

In this section, you will learn how to report your modelled effects in numbers. In general, this section will involve answering:

  • What are your modelled effects (with uncertainty)?

If you also have a categorical predictor with more than two levels, you will also want to answer:

  • Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)

Let’s examine each of these in turn:

5.2.1 What are your modelled effects (with uncertainty)?

Recall from the text above that your model effects describe the change in your response given a change in your predictor.

When you have a continuous predictor, this effect is the the amount of change in the response that is caused by a unit change in your continuous predictor. Information about this effect is captured in the coefficient (slope) of the linear fit. Expand the text box below to see why.

Here’s a model equation for a generalized linear model:

\[ \begin{align} g(\mu_i|Pred_i) &= \beta_1 \cdot Pred_i + \beta_0 \\ Resp_i &\sim F(\mu_i) \end{align} \]

where * \(g()\) is the link function * \(\mu_i\) is the mean fitted value \(\mu_i\) dependent on \(Pred_i\) * \(\beta_1\) is the coefficient of \(Pred_i\) * \(Pred_i\) is the value of your predictor * \(\beta_0\) is the intercept * \(Resp_i\) is the value of the response variable * \(F\) represents some error distribution assumption

As above, the effect of a continuous predictor on your response is the change in your fitted value (\(g(\mu_i)\)) for a unit change in your predictor (\(Pred_i\)):

\[ \begin{align} g(\mu|Pred+1) - g(\mu|Pred) &= (\beta_1 \cdot (Pred+1) + \beta_0) - (\beta_1 \cdot Pred_i + \beta_0) \\ &=\beta_1 \cdot Pred + \beta_1 + \beta_0 - \beta_1 \cdot Pred_i - \beta_0 \\ &=\beta_1 \end{align} \] For linear models with a normal error distribution assumption (i.e. when you are using link = "identity"), the link function \(g()\) is not there. For this reason, you can read the effects of your model directly from the model coefficients (\(\beta_1\), \(\beta_0\)) when you have a normal error distribution assumption. And for this reason, you will need to convert your model coefficients to quantify your modelled effects when you have an error distribution assumption that is not normal/uses a link function that isn’t “identity”.

This is because of the difference between the “link” scale and “response” scale for generalized linear models.

The link scale is where the model fitting happens, and where the evidence for your modelled effects is gathered.

The response scale is back in the real world - these are the actual effects in your response you would expect to see with a change in your predictor.

Recall that you can see the default link functions associated with each error distribution function with ?family:

When you have a categorical predictor, the coefficient represents the intercept of the linear fit and how that intercept changes as you move from one category to another. So if you have a categorical predictor with 3 categories (levels), you will have three coefficients that describe how the mean of your response changes as you move from category to category.

When you have an interaction among your predictors, you will have one coefficient for each combination in your interaction.

Finally, you need to report the uncertainty around your modelled effects so your peers know how much evidence there is for these modelled effects. There are a number of ways to do this. Two that we will use are i) reporting standard errors (SE) which are a measure of uncertainty of the average modelled effect, and ii) the 95% confidence intervals around your modelled effects, which tell you the range your modelled effects would be expected to fall into (with a 95% probability) if you redid your experiment.

These principles are true regardless of your model structure BUT our interpretation of what the effects mean in the real world will vary based on e.g. your error distribution assumptions (see the text box above for an explanation). In the examples below, you will learn how to quantify your modelled effects and uncertainty estimates for models with a normal error distribution assumption. There are sections below to show you how to do the same for models with other error distribution assumptions.

For a quick summary, here is how the method depends on your modelled error distribution assumption:

5.2.2 Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors)

When you have a categorical predictor in your best-specified model, you know there is evidence of an effect for at least one of your levels (categories). This step will tell you where you have evidence that the effects among your factor levels differ. The processes outlined below are the same for any model structure.

5.2.3 Models with a normal error distribution assumption:

Example 1: Resp1 ~ Cat1 + 1

Recall that Example 1 contains only one predictor and it is categorical:

Resp1 ~ Cat1 + 1

What are your modelled effects (with uncertainty)?

In the chapter on your Starting Model, you found your coefficients in the “Estimate” column of the summary() output of your model:

coef(summary(bestMod1)) # extract the coefficients from summary()
             Estimate Std. Error   t value     Pr(>|t|)
(Intercept) 20.954839   1.316231 15.920341 8.313153e-29
Cat1Sp2      5.876740   1.773637  3.313384 1.296521e-03
Cat1Sp3     -4.509677   1.861431 -2.422694 1.726174e-02

Interpreting the coefficients from this output takes practice - especially for categorical predictors because of the way that R treats categorical predictors in regression. With R’s “dummy coding”, one level of the predictor (here Sp1) is incorporated into the intercept estimate (21) as the reference level. The other coefficients in the Estimate column represent the change in modelled response when you move from the reference level (Sp1) to another level.

The modelled response when Cat1 = Sp2 is (5.9) units higher than this reference level = 21 + 5.9 = 11.8 units.

The modelled Resp1 when Cat1 = Sp3 is -4.5 units lower than the reference level = 21 + 5.9 = 1.4 units.5

So all the information we need is in this summary() output, but not easy to see immediately. An easier way is to use the emmeans6 package which helps you by reporting the coefficients directly. In our case, you use the emmeans package to get the mean value of the response at each level of the predictor. For categorical predictors, you do this with the emmeans() function:

library(emmeans) # load the emmeans package

emmOut <- emmeans(object = bestMod1, # your model
            specs = ~ Cat1, # your predictors
            type = "response") # report coefficients on the response scale

emmOut
 Cat1 emmean   SE df lower.CL upper.CL
 Sp1    21.0 1.32 97     18.3     23.6
 Sp2    26.8 1.19 97     24.5     29.2
 Sp3    16.4 1.32 97     13.8     19.1

Confidence level used: 0.95 

So now you have a modelled value of your response for each level of your categorical predictor - this captures the effect of your categorical predictor on your response. You also need to report uncertainty around this effect, and you can do this by reporting the standard error (SE) also reported by the emmeans() function.

When Cat1 is Sp1, Resp1 is estimated to be 21 ± 1.3 units.
When Cat1 is Sp2, Resp1 is estimated to be 26.8 ± 1.2 units.
When Cat1 is Sp3, Resp1 is estimated to be 16.4 ± 1.3 units.

Note that this is the same information in the summary() output just easier to read.7

Note also that two types of uncertainty are measured here. SE stands for the standard error around the prediction, and is a measure of uncertainty of the average modelled effect. The lower.CL and upper.CL represent the 95% confidence limits of the prediction - so if I observed a new Resp1 at a particular Cat1, there would be a 95% chance it would lie between the bounds of the confidence limits.

Finally, note that you can also get a quick plot of the effects by handing the emmeans() output to plot().

Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)

Since you have a categorical predictor with more than two levels, you will ask if all levels of Cat1 affect Resp1 in the same way - i.e. are the coefficients across factor levels significantly different from one another?

To get evidence about how each level affects your response, you need to test which effects differ from one another among the categories (levels) of your categorical predictor. This is done using multiple comparisons, i.e. you will compare the modelled effect of each level of your categorical predictor vs. each other level of your categorical predictor to determine which are different from each other (called pairwise testing).

This is done by testing the null hypothesis that the modelled effects of each pair are similar to one another, typically rejecting the null hypothesis when P < 0.05. Remember that the P-value is the probability that you observe a value at least as big as the one you observed even though our null hypothesis is true. In this case, you are looking at the value you are observing is the difference between coefficients estimated for two levels of your predictor. The P-value is the probability of getting a difference at least as big as the one you observed even though there is actually no difference between the coefficients (the null hypothesis is true).

A couple of things to note about multiple comparison testing:

  1. Multiple comparison testing on a categorical predictor should only be done after your hypothesis testing has given you evidence that the predictor has an effect on your response. That is, you should only do a multiple comparison test on a predictor if that predictor is in your best-specified model. As this is a test done after your hypothesis testing, it is called a post-hoc8 test.

  2. Multiple comparison testing can be a problem because you are essentially repeating a hypothesis test many times on the same data (i.e. are the effects of Sp1 different than those of Sp2? are the effects of Sp1 different than those of Sp3? are the effects of Sp2 different than those of Sp3?…). These repeated tests mean there is a high chance that you will find a difference in one test purely due to random chance, vs. due to there being an actual difference. To account for this, the multiple comparison tests you will perform have been formulated to correct for this increased error.

Multiple comparison testing is very simple with the emmeans package. It just requires you to hand the output from emmeans() to a function called pairs():

forComp <- pairs(emmOut)

forComp
 contrast  estimate   SE df t.ratio p.value
 Sp1 - Sp2    -5.88 1.77 97  -3.313  0.0037
 Sp1 - Sp3     4.51 1.86 97   2.423  0.0451
 Sp2 - Sp3    10.39 1.77 97   5.856  <.0001

P value adjustment: tukey method for comparing a family of 3 estimates 

The output shows the results of the multiple comparison (pairwise) testing. The values in the p.value column tell you the results of the hypothesis test comparing the coefficients between the two levels of your categorical predictor. For example, for Sp1 vs. Sp2, there is a 0.37% (p = 0.0037) probability of getting a difference in coefficients at least as big as 26.8 - 21 = r 5.8 even though the null hypothesis (no difference) is true. This value 0.37% (P = 0.0037) is small enough that you can say you have evidence that the effect of Cat1 on Resp1 is different for these two different levels (Sp1andSp2`).

Based on a threshold p-value of 0.05, we can see that:

There is evidence that the value of Resp1 when Cat1 is Sp1 is different (lower) than that when Cat1 is Sp2 as P = 0.0037 is less than P = 0.05.
There is a little evidence that the value of Resp1 when Cat1 is Sp1 is different (higher) than that when Cat1 is Sp3 as P = 0.0451 is less than P = 0.05 (but it is close!). There is some evidence that the value of Resp1 when Cat1 is Sp2 is different (higher) than that when Cat1 is Sp3 as P < 0.00019 is smaller than P = 0.05.10.

Note that the p-values have been adjusted via the Tukey method which adjusts the difference that the two coefficients need to have to allow for the fact that we are making multiple comparisons.11

Note that you can also get the results from the pairwise testing visually by handing the output of pairs() to plot().

Note that the difference between Resp1 when Cat1 is Sp1 vs. Sp3 is so close to overlapping zero. This is reflected in the high P value.

Example 2: Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1

Recall that Example 2 contains two predictors and both are categorical:

Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1

What are your modelled effects (with uncertainty)?

Again, for categorical predictors, there is a coefficient estimated for each level of the predictor. If there is one or more interactions among predictors, there will be one coefficient for each combination of levels among predictors. Let’s look at the summary() output of your model to understand better:

coef(summary(bestMod2)) # extract the coefficients from summary()
                        Estimate Std. Error    t value     Pr(>|t|)
(Intercept)            364.90879   24.73850 14.7506420 1.629062e-25
Cat2TypeB               37.37325   32.98467  1.1330491 2.602712e-01
Cat2TypeC              -75.94350   33.87459 -2.2419017 2.748208e-02
Cat2TypeD             -141.96558   38.32472 -3.7042819 3.695537e-04
Cat3Treat2             -43.62974   36.41409 -1.1981554 2.340735e-01
Cat3Control             64.11573   30.29835  2.1161457 3.715672e-02
Cat2TypeB:Cat3Treat2    42.70879   53.60009  0.7968045 4.277091e-01
Cat2TypeC:Cat3Treat2    73.60837   54.15228  1.3592849 1.775295e-01
Cat2TypeD:Cat3Treat2    66.02110   55.13228  1.1975037 2.343260e-01
Cat2TypeB:Cat3Control   81.25030   43.24327  1.8789121 6.356678e-02
Cat2TypeC:Cat3Control  -19.07730   41.29748 -0.4619483 6.452584e-01
Cat2TypeD:Cat3Control  -71.72444   46.17117 -1.5534463 1.239057e-01

Comparing this output to that of Example 1 shows many more estimated coefficients for Example 2. This is because we have one coefficient for each level of each predictor, as well as coefficients for each combination (the interaction term) of levels of the predictors.

Again, it takes a little practice to read the coefficients from the summary() output. For Example 2:

The modelled prediction for Resp2 when Cat2 is TypeA and Cat3 is Treat1 is 365 units (the intercept).

The modelled prediction for Resp2 when Cat2 is TypeB and Cat3 is Treat1 is 365 + 37 = 402 units.

The modelled prediction for Resp2 when Cat2 is TypeA and Cat3 is Treat2 is 365 - 44 = 321 units.

The modelled prediction for Resp2 when Cat2 is TypeC and Cat3 is Treat2 is 365 - 76 - 44 + 74 = 319 units.

As above, we can use the emmeans package to more easily see these coefficients:

emmOut <- emmeans(object = bestMod2, # your model
            specs = ~ Cat2 + Cat3 + Cat2:Cat3, # your predictors
            type = "response") # report coefficients on the response scale

emmOut
 Cat2  Cat3    emmean   SE df lower.CL upper.CL
 TypeA Treat1     365 24.7 88      316      414
 TypeB Treat1     402 21.8 88      359      446
 TypeC Treat1     289 23.1 88      243      335
 TypeD Treat1     223 29.3 88      165      281
 TypeA Treat2     321 26.7 88      268      374
 TypeB Treat2     401 32.7 88      336      466
 TypeC Treat2     319 32.7 88      254      384
 TypeD Treat2     245 29.3 88      187      304
 TypeA Control    429 17.5 88      394      464
 TypeB Control    548 21.8 88      504      591
 TypeC Control    334 15.9 88      302      366
 TypeD Control    215 18.9 88      178      253

Confidence level used: 0.95 

So now we have a modelled value of our response for each level of our categorical predictor. For example:

The modelled prediction for Resp2 when Cat2 is TypeA and Cat3 is Treat1 is 365 +/- 25 units.

The modelled prediction for Resp2 when Cat2 is TypeB and Cat3 is Treat1 is 402 +/- 22 units.

The modelled prediction for Resp2 when Cat2 is TypeA and Cat3 is Treat2 is 321 +/- 27 units.

The modelled prediction for Resp2 when Cat2 is TypeC and Cat3 is Treat2 is 319 +/- 33 units.

Finally, note that you can also get a quick plot of the effects by handing the emmeans() output to plot().

Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)

As with Example 1, you can also find out which combinations of predictor levels are leading to statistically different model predictions in Resp2:

forComp <- pairs(emmOut)

forComp
 contrast                      estimate   SE df t.ratio p.value
 TypeA Treat1 - TypeB Treat1    -37.373 33.0 88  -1.133  0.9923
 TypeA Treat1 - TypeC Treat1     75.944 33.9 88   2.242  0.5248
 TypeA Treat1 - TypeD Treat1    141.966 38.3 88   3.704  0.0180
 TypeA Treat1 - TypeA Treat2     43.630 36.4 88   1.198  0.9878
 TypeA Treat1 - TypeB Treat2    -36.452 41.0 88  -0.889  0.9991
 TypeA Treat1 - TypeC Treat2     45.965 41.0 88   1.120  0.9929
 TypeA Treat1 - TypeD Treat2    119.574 38.3 88   3.120  0.0940
 TypeA Treat1 - TypeA Control   -64.116 30.3 88  -2.116  0.6129
 TypeA Treat1 - TypeB Control  -182.739 33.0 88  -5.540  <.0001
 TypeA Treat1 - TypeC Control    30.905 29.4 88   1.051  0.9959
 TypeA Treat1 - TypeD Control   149.574 31.1 88   4.805  0.0004
 TypeB Treat1 - TypeC Treat1    113.317 31.8 88   3.563  0.0277
 TypeB Treat1 - TypeD Treat1    179.339 36.5 88   4.912  0.0002
 TypeB Treat1 - TypeA Treat2     81.003 34.5 88   2.348  0.4517
 TypeB Treat1 - TypeB Treat2      0.921 39.3 88   0.023  1.0000
 TypeB Treat1 - TypeC Treat2     83.338 39.3 88   2.119  0.6110
 TypeB Treat1 - TypeD Treat2    156.947 36.5 88   4.299  0.0025
 TypeB Treat1 - TypeA Control   -26.742 28.0 88  -0.956  0.9982
 TypeB Treat1 - TypeB Control  -145.366 30.9 88  -4.711  0.0005
 TypeB Treat1 - TypeC Control    68.278 27.0 88   2.531  0.3355
 TypeB Treat1 - TypeD Control   186.947 28.9 88   6.477  <.0001
 TypeC Treat1 - TypeD Treat1     66.022 37.3 88   1.769  0.8297
 TypeC Treat1 - TypeA Treat2    -32.314 35.3 88  -0.914  0.9988
 TypeC Treat1 - TypeB Treat2   -112.396 40.1 88  -2.804  0.1965
 TypeC Treat1 - TypeC Treat2    -29.979 40.1 88  -0.748  0.9998
 TypeC Treat1 - TypeD Treat2     43.631 37.3 88   1.169  0.9900
 TypeC Treat1 - TypeA Control  -140.059 29.0 88  -4.828  0.0003
 TypeC Treat1 - TypeB Control  -258.683 31.8 88  -8.134  <.0001
 TypeC Treat1 - TypeC Control   -45.038 28.1 88  -1.605  0.9027
 TypeC Treat1 - TypeD Control    73.631 29.9 88   2.465  0.3757
 TypeD Treat1 - TypeA Treat2    -98.336 39.6 88  -2.481  0.3654
 TypeD Treat1 - TypeB Treat2   -178.418 43.9 88  -4.064  0.0056
 TypeD Treat1 - TypeC Treat2    -96.001 43.9 88  -2.186  0.5636
 TypeD Treat1 - TypeD Treat2    -22.391 41.4 88  -0.541  1.0000
 TypeD Treat1 - TypeA Control  -206.081 34.1 88  -6.043  <.0001
 TypeD Treat1 - TypeB Control  -324.705 36.5 88  -8.894  <.0001
 TypeD Treat1 - TypeC Control  -111.061 33.3 88  -3.335  0.0532
 TypeD Treat1 - TypeD Control     7.609 34.8 88   0.218  1.0000
 TypeA Treat2 - TypeB Treat2    -80.082 42.2 88  -1.895  0.7587
 TypeA Treat2 - TypeC Treat2      2.335 42.2 88   0.055  1.0000
 TypeA Treat2 - TypeD Treat2     75.945 39.6 88   1.916  0.7460
 TypeA Treat2 - TypeA Control  -107.746 31.9 88  -3.374  0.0478
 TypeA Treat2 - TypeB Control  -226.369 34.5 88  -6.562  <.0001
 TypeA Treat2 - TypeC Control   -12.725 31.1 88  -0.409  1.0000
 TypeA Treat2 - TypeD Control   105.945 32.7 88   3.237  0.0693
 TypeB Treat2 - TypeC Treat2     82.417 46.3 88   1.781  0.8238
 TypeB Treat2 - TypeD Treat2    156.026 43.9 88   3.554  0.0285
 TypeB Treat2 - TypeA Control   -27.663 37.1 88  -0.745  0.9998
 TypeB Treat2 - TypeB Control  -146.287 39.3 88  -3.719  0.0172
 TypeB Treat2 - TypeC Control    67.357 36.4 88   1.852  0.7846
 TypeB Treat2 - TypeD Control   186.027 37.8 88   4.923  0.0002
 TypeC Treat2 - TypeD Treat2     73.609 43.9 88   1.677  0.8739
 TypeC Treat2 - TypeA Control  -110.081 37.1 88  -2.967  0.1366
 TypeC Treat2 - TypeB Control  -228.704 39.3 88  -5.815  <.0001
 TypeC Treat2 - TypeC Control   -15.060 36.4 88  -0.414  1.0000
 TypeC Treat2 - TypeD Control   103.609 37.8 88   2.742  0.2240
 TypeD Treat2 - TypeA Control  -183.690 34.1 88  -5.387  <.0001
 TypeD Treat2 - TypeB Control  -302.313 36.5 88  -8.281  <.0001
 TypeD Treat2 - TypeC Control   -88.669 33.3 88  -2.663  0.2624
 TypeD Treat2 - TypeD Control    30.000 34.8 88   0.861  0.9993
 TypeA Control - TypeB Control -118.624 28.0 88  -4.242  0.0030
 TypeA Control - TypeC Control   95.021 23.6 88   4.023  0.0064
 TypeA Control - TypeD Control  213.690 25.7 88   8.299  <.0001
 TypeB Control - TypeC Control  213.644 27.0 88   7.918  <.0001
 TypeB Control - TypeD Control  332.314 28.9 88  11.514  <.0001
 TypeC Control - TypeD Control  118.669 24.7 88   4.809  0.0004

P value adjustment: tukey method for comparing a family of 12 estimates 

This is a hard table to navigate as every possible combination of the levels across predictors is compared. This may not be what you want. Instead, you can look for differences of one predictor based on a value of another. For example:

forComp <- pairs(emmOut, # emmeans output
                 simple = "Cat2") # contrasts within this categorical predictor

forComp
Cat3 = Treat1:
 contrast      estimate   SE df t.ratio p.value
 TypeA - TypeB   -37.37 33.0 88  -1.133  0.6702
 TypeA - TypeC    75.94 33.9 88   2.242  0.1201
 TypeA - TypeD   141.97 38.3 88   3.704  0.0021
 TypeB - TypeC   113.32 31.8 88   3.563  0.0033
 TypeB - TypeD   179.34 36.5 88   4.912  <.0001
 TypeC - TypeD    66.02 37.3 88   1.769  0.2949

Cat3 = Treat2:
 contrast      estimate   SE df t.ratio p.value
 TypeA - TypeB   -80.08 42.2 88  -1.895  0.2375
 TypeA - TypeC     2.34 42.2 88   0.055  0.9999
 TypeA - TypeD    75.94 39.6 88   1.916  0.2288
 TypeB - TypeC    82.42 46.3 88   1.781  0.2894
 TypeB - TypeD   156.03 43.9 88   3.554  0.0034
 TypeC - TypeD    73.61 43.9 88   1.677  0.3422

Cat3 = Control:
 contrast      estimate   SE df t.ratio p.value
 TypeA - TypeB  -118.62 28.0 88  -4.242  0.0003
 TypeA - TypeC    95.02 23.6 88   4.023  0.0007
 TypeA - TypeD   213.69 25.7 88   8.299  <.0001
 TypeB - TypeC   213.64 27.0 88   7.918  <.0001
 TypeB - TypeD   332.31 28.9 88  11.514  <.0001
 TypeC - TypeD   118.67 24.7 88   4.809  <.0001

P value adjustment: tukey method for comparing a family of 4 estimates 

Here you can more easily see the contrasts, and the adjustment to the P-value will be limited to what you actually need. Based on a threshold P-value of 0.05, you can see that the effects at some combinations of your predictors are not statistically different from each other.

For example, the coefficient (predicted Resp2) when Cat3 = Treat1 and Cat2 = TypeA is 365 and this is 37 lower than the coefficient when Cat3 = Treat1 and Cat2 = TypeB, but there is no evidence that the difference has any meaning as P = 0.67 for this comparison.

On the other hand, some combinations of your predictors are statistically different from each other. For example, a comparison of modelled Resp2 when Cat3 = Control and Cat2 = TypeB (548) vs. Cat3 = Control and Cat2 = TypeD (215) shows that they differ by 333 and that there is strong evidence that this difference is meaningful as P < 0.001.

Note that you could also organize this with contrasts by Cat3 instead:

forComp <- pairs(emmOut, # emmeans output
                 simple = "Cat3") # contrasts within this categorical predictor

forComp
Cat2 = TypeA:
 contrast         estimate   SE df t.ratio p.value
 Treat1 - Treat2    43.630 36.4 88   1.198  0.4574
 Treat1 - Control  -64.116 30.3 88  -2.116  0.0924
 Treat2 - Control -107.745 31.9 88  -3.374  0.0031

Cat2 = TypeB:
 contrast         estimate   SE df t.ratio p.value
 Treat1 - Treat2     0.921 39.3 88   0.023  0.9997
 Treat1 - Control -145.366 30.9 88  -4.711  <.0001
 Treat2 - Control -146.287 39.3 88  -3.719  0.0010

Cat2 = TypeC:
 contrast         estimate   SE df t.ratio p.value
 Treat1 - Treat2   -29.979 40.1 88  -0.748  0.7357
 Treat1 - Control  -45.038 28.1 88  -1.605  0.2489
 Treat2 - Control  -15.060 36.4 88  -0.414  0.9099

Cat2 = TypeD:
 contrast         estimate   SE df t.ratio p.value
 Treat1 - Treat2   -22.391 41.4 88  -0.541  0.8514
 Treat1 - Control    7.609 34.8 88   0.218  0.9741
 Treat2 - Control   30.000 34.8 88   0.861  0.6661

P value adjustment: tukey method for comparing a family of 3 estimates 

Note that you can also get the results from the pairwise testing visually by handing the output of pairs() to plot().

Example 3: Resp3 ~ Cont4 + 1

Recall that Example 3 contains one predictor and the predictor is continuous:

Resp3 ~ Cont4 + 1

What are your modelled effects (with uncertainty)?

For a continuous predictor, the effect is captured in one coefficient that describes the change in the response for a unit change in the continuous predictor. For models with a normal error distribution assumption, this is communicated as the slope.

Let’s look at the summary() output for Example 3:

coef(summary(bestMod3)) # extract the coefficients from summary()
             Estimate Std. Error   t value     Pr(>|t|)
(Intercept) -226.9131  166.09360 -1.366176 1.750105e-01
Cont4        260.0574   14.56593 17.853818 1.438227e-32

This tells you that for a unit change in Cont4, you get a 260 ± 1512.

Note that this same information is in the emmeans package. Because Cont4 is a continuous predictor, you need the emtrends() function, which provides the 95% confidence interval on the estimate:

trendsOut <- emtrends(bestMod3,  # your model
                      specs = ~ Cont4 + 1, # your predictors
                      var = "Cont4") # your continuous predictors

trendsOut
 Cont4 Cont4.trend   SE df lower.CL upper.CL
  9.94         260 14.6 98      231      289

Confidence level used: 0.95 

Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)

This question is only applies to categorical predictors of which there are none in Example 3.

Example 4: Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1

Recall that Example 4 contains two predictors, one is categorical and one is continuous:

Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1

What are your modelled effects (with uncertainty)?

The effects estimated for this model will include values of the response (Resp4) at each level of the categorical predictor (Cat5) as well as slopes describing the change in the response when the continuous predictor (Cont6) changes, and a different slope will be estimated for each level in Cat5 (this is because of the interaction in the model).

As above, you can take a look at your modelled coefficients with the summary() output:

coef(summary(bestMod4))
                   Estimate  Std. Error     t value     Pr(>|t|)
(Intercept)     98.34679337 1.868900753 52.62280151 1.536692e-71
Cat5Urban       -0.23642403 2.895404981 -0.08165491 9.350947e-01
Cat5Wild        -0.86733554 3.238562391 -0.26781498 7.894285e-01
Cont6           -0.00293078 0.003740830 -0.78345704 4.353283e-01
Cat5Urban:Cont6  0.01511745 0.005611511  2.69400763 8.361918e-03
Cat5Wild:Cont6   0.03085821 0.006183238  4.99062262 2.756559e-06

This shows:

The model prediction of Resp4 when Cat5 is Farm and Cont6 is 0 is 98.34 units13 (the intercept).

The model prediction of Resp4 when Cat5 is Urban and Cont6 is 0 is 98.34 - 0.24 = 98.1 units.

The model prediction of Resp4 when Cat5 is Wild and Cont6 is 0 is 98.34 - 0.87 = 97.5 units.

The slope of the relationship between Cont6 and Resp4 when Cat5 is Farm is -0.0029.14

The slope of the relationship between Cont6 and Resp4 when Cat5 is Urban is -0.0029 + 0.015 = 0.0121.

The slope of the relationship between Cont6 and Resp4 when Cat5 is Wild is -0.0029 + 0.031 = 0.0281.

Again, interpreting the coefficients from the summary() output is tedious and not necessary: You can use the emmeans package to give you the modelled response for each level of the categorical predictor (Cat5) directly:

emmOut <- emmeans(object = bestMod4, # your model
            specs = ~ Cat5 + Cont6 + Cat5:Cont6 + 1, # your predictors
            type = "response") # report coefficients on the response scale

emmOut
 Cat5  Cont6 emmean    SE df lower.CL upper.CL
 Farm    510   96.9 0.458 94     95.9     97.8
 Urban   510  104.3 0.421 94    103.5    105.2
 Wild    510  111.7 0.511 94    110.7    112.7

Confidence level used: 0.95 

Note that emmeans() sets our continuous predictor (Cont6) to the mean value of Cont6 (510 units). We can also control this with the at = argument:

emmOut <- emmeans(object = bestMod4, # your model
            specs = ~ Cat5 + Cont6 + Cat5:Cont6 + 1, # your predictors
            type = "response", # report coefficients on the response scale
            at = list(Cont6 = 0)) # control the value of your continuous predictor at which to make the coefficient estimates


emmOut
 Cat5  Cont6 emmean   SE df lower.CL upper.CL
 Farm      0   98.3 1.87 94     94.6      102
 Urban     0   98.1 2.21 94     93.7      103
 Wild      0   97.5 2.64 94     92.2      103

Confidence level used: 0.95 

By setting at = 0, you get the intercept - i.e. the modelled Resp4 when Cont6 = 0 for each level of Cat5, and this is what is reported in the summary() output.

Similarly, you can get the emmeans package to give you the slope coefficients for the continuous predictor (Cont6) using the emtrends() function, rather than interpreting them from the summary() output:

trendsOut <- emtrends(bestMod4,  # your model
                      specs = ~ Cat5 + Cont6 + Cat5:Cont6 + 1, # your predictors
                      var = "Cont6") # your continuous predictors


trendsOut
 Cat5  Cont6 Cont6.trend      SE df lower.CL upper.CL
 Farm    510    -0.00293 0.00374 94 -0.01036   0.0045
 Urban   510     0.01219 0.00418 94  0.00388   0.0205
 Wild    510     0.02793 0.00492 94  0.01815   0.0377

Confidence level used: 0.95 

Note that the 95% confidence limits for the modelled effect for Cont6 when Cat5 = Farm includes zero. We can also see this with the test() function.

test(trendsOut) # get test if coefficient is different than zero.
 Cat5  Cont6 Cont6.trend      SE df t.ratio p.value
 Farm    510    -0.00293 0.00374 94  -0.783  0.4353
 Urban   510     0.01219 0.00418 94   2.914  0.0045
 Wild    510     0.02793 0.00492 94   5.673  <.0001

combining the two output, we can report that:

  • there is no evidence on an effect on Resp4 of Cont6 when Cat5 = Farm (p = 0.44).

  • there is evidence for a positive effect of Cont6 on Resp4 when Cat5 = Urban. The effect is 0.01215 with a 95% confidence interval of 0.0039 - 0.020516.

  • there is evidence for a positive effect of Cont6 on Resp4 when Cat5 = Wild. The effect is 0.02817 with a 95% confidence interval of 0.0182 - 0.037718.

Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)

Since you have a categorical predictor with more than two levels, you will want to report where the effects among levels differ from one another.

Note that the effects in Example 4 include slope coefficients (associated with the continuous predictor) and intercept coefficients (associated with the categorical predictor). Because the slope and intercept are dependent on one another (when the slope changes, the intercept will naturally change), you need to first check if the slope coefficients vary across your categorical levels. Only if they don’t vary (i.e. the modelled effects produce lines that are parallel across your categorical levels) would you go on to test if the intercept coefficiens vary across your categorical levels.

You can find out which slopes (i.e. the effect of Cont6 on Resp4) are different across the levels of Cat5 using emtrends() with a request for pairwise testing:

forCompSlope <- pairs(trendsOut)

forCompSlope
 contrast                               estimate      SE df t.ratio p.value
 Farm Cont6509.767 - Urban Cont6509.767  -0.0151 0.00561 94  -2.694  0.0226
 Farm Cont6509.767 - Wild Cont6509.767   -0.0309 0.00618 94  -4.991  <.0001
 Urban Cont6509.767 - Wild Cont6509.767  -0.0157 0.00646 94  -2.437  0.0437

P value adjustment: tukey method for comparing a family of 3 estimates 

Here you see that the slope estimates (the effect of Cont6 on Resp4) for all levels of Cat5 are significantly different from one another (0.0001 ≤ P ≤ 0.044).

Note that you can also get the results from the pairwise testing visually by handing the output of pairs() to plot().

Again: only if all the slopes were similar, would you want to test if the levels of your categorical predictor (Cat4) have significantly different coefficients (intercepts) from each other. You would do this with pairs() on the output from emmeans() (the emmOut object).


5.2.4 Models with a poisson error distribution assumption:

Background

For models with a poisson error distribution assumption, you will typically start with a link function that is the natural logarithm or link = "log". You need to take this link function into consideration when you report your modelled effects.

If you are using link = "log" (as with a poisson error distribution, and sometimes also a Gamma error distribution - more on this to come), you get your coefficient on the response scale by taking e to the coefficient on the link scale (with the R function exp()). This coefficient is called the rate ratio19 and it tells you the % change in the response for a unit change in your predictor.


For a continuous predictor, the effect is captured in one coefficient that describes the change in the response for a unit change in the continuous predictor.

For models with a normal error distribution assumption, this was simply a slope.

For models using a log link (including many models with a poisson error distribution assumption), this is communicated as the rate ratio or the % change in the response for a unit change in your predictor. This allows you to communicate the effect of the continuous predictor while accounting for the curve in the relationship (see visualization of model effects above). This curve occurs because the model is not allowed to go below zero to be consistent with a poisson error distribution assumption.

But why do we convert coefficients from the link to response scale with \(e^x\)20 when link = "log"?

A reminder that we can present our linear model like this:

\[ \begin{align} g(\mu_i|Pred_i) &= \beta_1 \cdot Pred_i + \beta_0 \\ Resp_i &\sim F(\mu_i) \end{align} \]

where * \(g()\) is the link function * \(\mu_i\) is the mean fitted value \(\mu_i\) dependent on \(Pred_i\) * \(\beta_1\) is the coefficient of \(Pred_i\) * \(Pred_i\) is the value of your predictor * \(\beta_0\) is the intercept * \(Resp_i\) is the value of the response variable * \(F\) represents some error distribution assumption

If you have an error distribution assumption (\(F\)), the link function \(g()\) is \(log_e\)21:

\[ \begin{align} log_e(\mu_i|Pred_i) &= \beta_1 \cdot Pred_i + \beta_0 \\ Resp_i &\sim poisson(\mu_i) \end{align} \]

Then the rate ratio (RR), or % change in the response for a unit change in the predictor becomes,

\[ \begin{align} RR&=\frac{(\mu_i|Pred = a+1)}{(\mu_i|Pred = a)}\\[2em] log_e(RR)&=log_e(\frac{(\mu_i|Pred = a+1)}{(\mu_i|Pred = a)})\\[2em] log_e(RR)&=log_e(\mu_i|Pred = a+1)-log_e(\mu_i|Pred = a)\\[2em] log_e(RR)&=(\beta_1\cdot (a+1)+\beta_0)-(\beta_1\cdot (a)+\beta_0)\\[2em] log_e(RR)&=\beta_1\cdot a + \beta_1+\beta_0-\beta_1\cdot a-\beta_0\\[2em] log_e(RR)&=\beta_1\\[2em] RR &=e^{\beta_1} \end{align} \]

One thing to note:

The emmeans() function will convert the coefficients (intercepts) from the link to the response scale. You can ask for this with the type = "response" argument.

In contrast, the emtrends() function does not convert the coefficients (slopes) to represent effects on the response scale. This is because emtrends() is reporting the slope of a straight line - the trend line on the link scale. But the line isn’t straight on the response scale.

We need to convert from the link to the response by hand when we use emtrends()^[remember, the coefficients for categorical predictors are also being converted from the link to the response scale. It is just that R is doing it for us in the emmeans() function when we add the argument `type = “response”). How we make the conversion, and interpret the result, will depend on our error distribution assumption:

Here we will go through examples of reporting for models with a poisson error distribution assumption. If you want more details on why you are using the methods below, check the section on “Models with a normal error distribution assumption” (Section 5.2.3) for context.

Examples

Here I have made new examples where the error distribution assumption is poisson. This changes the response variables, indicated by Resp1.p, Resp2.p, etc. to be integers ≥ 0.

Example 1: Resp1.p ~ Cat1 + 1

Recall that Example 1 contains only one predictor and it is categorical:

Resp1.p ~ Cat1 + 1

The best-specified model (now with poisson error distribution assumption) is:

bestMod1.p

Call:  glm(formula = Resp1.p ~ Cat1 + 1, family = poisson(link = "log"), 
    data = Dat)

Coefficients:
(Intercept)      Cat1StB      Cat1StC  
    6.15524     -0.04869      0.09948  

Degrees of Freedom: 49 Total (i.e. Null);  47 Residual
Null Deviance:      139 
Residual Deviance: 40.01    AIC: 446

that was fit to data in myDat1.p:

str(myDat1.p)
'data.frame':   50 obs. of  2 variables:
 $ Cat1   : Factor w/ 3 levels "StA","StB","StC": 2 2 2 1 3 2 3 1 3 1 ...
 $ Resp1.p: int  434 476 443 473 516 449 529 462 512 429 ...

and the dredge() table you used to pick your bestMod1.p is in dredgeOut1.p

dredgeOut1.p
Global model call: glm(formula = Resp1.p ~ Cat1 + 1, family = poisson(link = "log"), 
    data = Dat)
---
Model selection table 
  (Intrc) Cat1 df   logLik  AICc delta weight
2   6.155    +  3 -219.983 446.5  0.00      1
1   6.164       1 -269.477 541.0 94.55      0
Models ranked by AICc(x) 

And here’s a visualization of the model effects:

# Set up your predictors for the visualized fit
forCat1<-unique(myDat1.p$Cat1) # every value of your categorical predictor

# create a data frame with your predictors
forVis<-expand.grid(Cat1=forCat1) # expand.grid() function makes sure you have all combinations of predictors

# Get your model fit estimates at each value of your predictors
modFit<-predict(bestMod1.p, # the model
                newdata = forVis, # the predictor values
                type = "response", # make the predictions on the response scale
                se.fit = TRUE) # include uncertainty estimate

forVis$Fit<-modFit$fit # add your fit to the data frame
forVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frame
forVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame


library(ggplot2)

ggplot() + # start ggplot
  geom_point(data = myDat1.p, # add observations to your plot
             mapping = aes(x = Cat1, y = Resp1.p), 
             position=position_jitter(width=0.1)) + # control position of data points so they are easier to see on the plot
  geom_errorbar(data = forVis, # add the uncertainty to your plot
              mapping = aes(x = Cat1, y = Fit, ymin = Lower, ymax = Upper),
              linewidth=1.2) + # control thickness of errorbar line
  geom_point(data = forVis, # add the modelled fit to your plot
              mapping = aes(x = Cat1, y = Fit), 
             shape = 22, # shape of point
             size = 3, # size of point
             fill = "white", # fill colour for plot
             col = 'black') + # outline colour for plot
  ylab("Resp1.p, (units)")+ # change y-axis label
  xlab("Cat1, (units)")+ # change x-axis label
  theme_bw() # change ggplot theme

What are your modelled effects (with uncertainty)?

For categorical predictors, you can get the coefficients on the response scale directly with the emmeans() function22 when you specify type = “response”:

library(emmeans) # load the emmeans package

emmOut <- emmeans(object = bestMod1.p, # your model
            specs = ~ Cat1 + 1, # your predictors
            type = "response") # report coefficients on the response scale

emmOut
 Cat1 rate   SE  df asymp.LCL asymp.UCL
 StA   471 6.54 Inf       459       484
 StB   449 4.32 Inf       440       457
 StC   520 5.89 Inf       509       532

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

Notice that the column reported is “rate”: these are the coefficients converted back onto the response scale.

Note also that two types of uncertainty are measured here. SE stands for the standard error around the prediction, and is a measure of uncertainty of the average modelled effect. The lower.CL and upper.CL represent the 95% confidence limits of the prediction - so if I observed a new Resp1.p at a particular Cat1, there would be a 95% chance it would lie between the bounds of the confidence limits.

From this output you see:

When Cat1 is SpA, Resp1.p is estimated to be 471 (95% confidence interval: 459 to 484 units).
When Cat1 is SpB, Resp1.p is estimated to be 449 (95% confidence interval: 440 to 457 units). When Cat1 is SpC, Resp1.p is estimated to be 520 (95% confidence interval: 509 to 532 units).

Finally, note that you can also get a quick plot of the effects by handing the emmeans() output to plot().

Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)

Here you can compare the effects of the levels of Cat1 on Resp1.p to see if the effects differ between the levels.

forComp <- pairs(emmOut)

forComp
 contrast  ratio     SE  df null z.ratio p.value
 StA / StB 1.050 0.0177 Inf    1   2.880  0.0111
 StA / StC 0.905 0.0162 Inf    1  -5.552  <.0001
 StB / StC 0.862 0.0128 Inf    1  -9.968  <.0001

P value adjustment: tukey method for comparing a family of 3 estimates 
Tests are performed on the log scale 

The output shows the results of the multiple comparison (pairwise) testing.

The values in the p.value column tell you the results of the hypothesis test comparing the coefficients between the two levels of your categorical predictor.

For example:

  • there is evidence that Resp1.p is significantly larger when Cat 1 = SpA than when Cat1 = SpB (p = 0.0111). Resp1.p is expected to be ~ 5 ± 2% bigger when Cat 1 = SpA than when Cat 1 = SpB.

  • there is evidence that Resp1.p is significantly smaller when Cat 1 = SpA than when Cat1 = SpC (p = 0). Resp1.p is expected to be ~ 9 ± 2% smaller when Cat 1 = SpA than when Cat 1 = SpC.

Note that you can also get the results from the pairwise testing visually by handing the output of pairs() to plot().

Example 2: Resp2.p ~ Cat2 + Cat3 + Cat2:Cat3 + 1

Recall that Example 2 contains two predictors and both are categorical:

Resp2.p ~ Cat2 + Cat3 + Cat2:Cat3 + 1

The best-specified model (now with poisson error distribution assumption) is:

bestMod2.p

Call:  glm(formula = Resp2.p ~ Cat2 + Cat3 + Cat2:Cat3 + 1, family = poisson(link = "log"), 
    data = myDat2.p)

Coefficients:
        (Intercept)            Cat2TypeB            Cat2TypeC  
             2.7213              -0.3365              -0.4526  
          Cat2TypeD            Cat3Treat  Cat2TypeB:Cat3Treat  
            -1.2744              -0.3810              -0.6176  
Cat2TypeC:Cat3Treat  Cat2TypeD:Cat3Treat  
            -0.5999               0.4381  

Degrees of Freedom: 49 Total (i.e. Null);  42 Residual
Null Deviance:      132.5 
Residual Deviance: 35.27    AIC: 239.5

that was fit to data in myDat2.p:

str(myDat2.p)
'data.frame':   50 obs. of  3 variables:
 $ Cat2   : Factor w/ 4 levels "TypeA","TypeB",..: 1 2 4 1 2 4 1 1 2 2 ...
 $ Cat3   : Factor w/ 2 levels "Control","Treat": 2 2 1 1 2 1 1 1 1 1 ...
 $ Resp2.p: int  10 3 3 19 4 4 9 13 15 12 ...

and the dredge() table you used to pick your bestMod2.p is in

dredgeOut2.p
Global model call: glm(formula = Resp2.p ~ Cat2 + Cat3 + Cat2:Cat3 + 1, family = poisson(link = "log"), 
    data = myDat2.p)
---
Model selection table 
  (Int) Ct2 Ct3 Ct2:Ct3 df   logLik  AICc delta weight
8 2.721   +   +       +  8 -111.758 243.0  0.00   0.77
4 2.860   +   +          5 -117.038 245.4  2.41   0.23
2 2.461   +              4 -134.275 277.4 34.41   0.00
3 2.344       +          2 -150.075 304.4 61.38   0.00
1 2.087                  1 -160.376 322.8 79.81   0.00
Models ranked by AICc(x) 

And here’s a visualization of the model effects:

#### i) choosing the values of your predictors at which to make a prediction

# Set up your predictors for the visualized fit
forCat2<-unique(myDat2.p$Cat2) # every level of your categorical predictor
forCat3<-unique(myDat2.p$Cat3) # every level of your categorical predictor
  
# create a data frame with your predictors
forVis<-expand.grid(Cat2 = forCat2, Cat3 = forCat3) # expand.grid() function makes sure you have all combinations of predictors

#### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor


# Get your model fit estimates at each value of your predictors
modFit<-predict(bestMod2.p, # the model
                newdata = forVis, # the predictor values
                type = "response", # make the predictions on the response scale
                se.fit = TRUE) # include uncertainty estimate

forVis$Fit<-modFit$fit # add your fit to the data frame
forVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frame
forVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame

#### iii) use the model estimates to plot your model fit

library(ggplot2) # load ggplot2 library

ggplot() + # start ggplot
  geom_point(data = myDat2.p, # add observations to your plot
             mapping = aes(x = Cat2, y = Resp2.p, col = Cat3), 
             position=position_jitterdodge(jitter.width=0.75, dodge.width=0.75)) + # control position of data points so they are easier to see on the plot
  geom_errorbar(data = forVis, # add the uncertainty to your plot
              mapping = aes(x = Cat2, y = Fit, ymin = Lower, ymax = Upper, col = Cat3),
              position=position_dodge(width=0.75), # control position of data points so they are easier to see on the plot
              size=1.2) + # control thickness of errorbar line
  geom_point(data = forVis, # add the modelled fit to your plot
              mapping = aes(x = Cat2, y = Fit, fill = Cat3), 
             shape = 22, # shape of point
             size=3, # size of point
             col = 'black', # outline colour for point
             position=position_dodge(width=0.75)) + # control position of data points so they are easier to see on the plot
  ylab("Resp2.p, (units)")+ # change y-axis label
  xlab("Cat2, (units)")+ # change x-axis label
  labs(fill="Cat3, (units)", col="Cat3, (units)") + # change legend title
  theme_bw() # change ggplot theme

What are your modelled effects (with uncertainty)?

As above, we can use the emmeans package to see the effects of the predictors on the scale of the response:

emmOut <- emmeans(object = bestMod2.p, # your model
            specs = ~ Cat2 + Cat3 + Cat2:Cat3, # your predictors
            type = "response") # report coefficients on the response scale

emmOut
 Cat2  Cat3     rate    SE  df asymp.LCL asymp.UCL
 TypeA Control 15.20 1.740 Inf     12.14     19.03
 TypeB Control 10.86 1.250 Inf      8.67     13.59
 TypeC Control  9.67 1.800 Inf      6.72     13.91
 TypeD Control  4.25 1.030 Inf      2.64      6.84
 TypeA Treat   10.38 0.894 Inf      8.77     12.29
 TypeB Treat    4.00 0.707 Inf      2.83      5.66
 TypeC Treat    3.62 0.673 Inf      2.52      5.22
 TypeD Treat    4.50 1.500 Inf      2.34      8.65

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

So now we have a modelled value of our response for each level of our categorical predictors. For example:

The modelled prediction for Resp2.p when Cat2 is TypeA and Cat3 is Treat is 15.2 (95% confidence interval: 12.1 - 19.0 units).

The modelled prediction for Resp2.p when Cat2 is TypeC and Cat3 is Control is 9.7 units (95% confidence interval: 6.7 - 13.9 units).

Finally, note that you can also get a quick plot of the effects by handing the emmeans() output to plot().

Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)

As with Example 1, you can also find out which combinations of predictor levels are leading to statistically different model predictions in Resp2.p:

forComp <- pairs(emmOut, # emmeans output
                 simple = "Cat2") # contrasts within this categorical predictor

forComp
Cat3 = Control:
 contrast      ratio    SE  df null z.ratio p.value
 TypeA / TypeB 1.400 0.227 Inf    1   2.074  0.1616
 TypeA / TypeC 1.572 0.343 Inf    1   2.074  0.1617
 TypeA / TypeD 3.576 0.960 Inf    1   4.750  <.0001
 TypeB / TypeC 1.123 0.245 Inf    1   0.532  0.9513
 TypeB / TypeD 2.555 0.685 Inf    1   3.496  0.0027
 TypeC / TypeD 2.275 0.695 Inf    1   2.690  0.0359

Cat3 = Treat:
 contrast      ratio    SE  df null z.ratio p.value
 TypeA / TypeB 2.596 0.510 Inf    1   4.852  <.0001
 TypeA / TypeC 2.865 0.586 Inf    1   5.142  <.0001
 TypeA / TypeD 2.308 0.794 Inf    1   2.429  0.0717
 TypeB / TypeC 1.103 0.283 Inf    1   0.384  0.9808
 TypeB / TypeD 0.889 0.335 Inf    1  -0.312  0.9895
 TypeC / TypeD 0.806 0.307 Inf    1  -0.567  0.9420

P value adjustment: tukey method for comparing a family of 4 estimates 
Tests are performed on the log scale 

Here you can more easily see the contrasts. Contrast ratios similar to 1 mean that the predictor levels lead to a similar value of the response.

Based on a threshold P-value of 0.05, you can see that the effects at some combinations of your predictors are not statistically different from each other. For example, the coefficient (predicted Resp2.p) when Cat3 = Control and Cat2 = TypeA (15.2) is 1.4 times (or 40%) more than the predicted Resp2.p when Cat3 = Control and Cat2 = TypeB (10.9), but there is no evidence that the difference has any meaning as P = 0.16 for this comparison.

On the other hand, some combinations of your predictors are statistically different from each other. For example, a comparison of modelled Resp2.p when Cat3 = Control and Cat2 = TypeB (10.9) is 2.6 times (or 160%) higher than when Cat3 = Control and Cat2 = TypeD (4.2) and there is strong evidence that this difference is meaningful as P < 0.0027.

Note that you can also get the results from the pairwise testing visually by handing the output of pairs() to plot().

Example 3: Resp3.p ~ Cont4 + 1

Recall that Example 3 contains one predictor and the predictor is continuous:

Resp3 ~ Cont4 + 1

The best-specified model (now with poisson error distribution assumption) is:

bestMod3.p

Call:  glm(formula = Resp ~ Cont4 + 1, family = poisson(link = "log"), 
    data = myDat3.p)

Coefficients:
(Intercept)        Cont4  
   1.577204     0.006674  

Degrees of Freedom: 49 Total (i.e. Null);  48 Residual
Null Deviance:      74.74 
Residual Deviance: 62.96    AIC: 256.7

that was fit to data in myDat3.p:

str(myDat3.p)
'data.frame':   50 obs. of  2 variables:
 $ Cont4  : num  67.5 61.9 71.5 116.6 55.9 ...
 $ Resp3.p: int  12 8 4 12 0 2 4 10 8 7 ...

and the dredge() table you used to pick your bestMod3.p is

dredgeOut3.p
Global model call: glm(formula = Resp ~ Cont4 + 1, family = poisson(link = "log"), 
    data = myDat3.p)
---
Model selection table 
  (Intrc)    Cont4 df   logLik  AICc delta weight
2   1.577 0.006674  2 -126.339 256.9  0.00  0.992
1   2.069           1 -132.230 266.5  9.61  0.008
Models ranked by AICc(x) 

And here’s a visualization of the model effects:

#### i) choosing the values of your predictors at which to make a prediction


# Set up your predictors for the visualized fit
forCont4<-seq(from = min(myDat3.p$Cont4), to = max(myDat3.p$Cont4), by = 1)# a sequence of numbers in your continuous predictor range
  
# create a data frame with your predictors
forVis<-expand.grid(Cont4 = forCont4) # expand.grid() function makes sure you have all combinations of predictors.  

#### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor


# Get your model fit estimates at each value of your predictors
modFit<-predict(bestMod3.p, # the model
                newdata = forVis, # the predictor values
                type = "response", # make the predictions on the response scale
                se.fit = TRUE) # include uncertainty estimate

forVis$Fit<-modFit$fit # add your fit to the data frame
forVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frame
forVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame

#### iii) use the model estimates to plot your model fit


library(ggplot2) # load ggplot2 library

ggplot() + # start ggplot
  
  geom_point(data = myDat3.p, # add observations to your plot
             mapping = aes(x = Cont4, y = Resp3.p)) + # control position of data points so they are easier to see on the plot
  
  geom_line(data = forVis, # add the modelled fit to your plot
              mapping = aes(x = Cont4, y = Fit),
              size = 1.2) + # control thickness of line
  
    geom_line(data = forVis, # add uncertainty to your plot (upper line)
              mapping = aes(x = Cont4, y = Upper),
              size = 0.8, # control thickness of line
              linetype = 2) + # control style of line
  
      geom_line(data = forVis, # add uncertainty to your plot (lower line)
              mapping = aes(x = Cont4, y = Lower),
              size = 0.8, # control thickness of line
              linetype = 2) + # control style of line
  
  ylab("Resp3.p, (units)") + # change y-axis label
  
  xlab("Cont4, (units)") + # change x-axis label
  
  theme_bw() # change ggplot theme

What are your modelled effects (with uncertainty)?

trendsOut <- emtrends(bestMod3.p,
                      specs = ~ Cont4 + 1, # your predictors
                      var = "Cont4") # your predictor of interest

trendsOut
 Cont4 Cont4.trend      SE  df asymp.LCL asymp.UCL
  71.5     0.00667 0.00193 Inf   0.00288    0.0105

Confidence level used: 0.95 

Note that these are the same estimates as we find in the summary() of your model object

coef(summary(bestMod3.p))
               Estimate  Std. Error   z value     Pr(>|z|)
(Intercept) 1.577203770 0.155306248 10.155443 3.133832e-24
Cont4       0.006674037 0.001934204  3.450534 5.594789e-04

This is an estimate of the modelled effects in the linked world. We need to consider our log link function to report the modelled effects in the response world.

Since you have a poisson error distribution assumption, you can convert the estimate made by emtrends() on to the response scale with:

trendsOut.df <- data.frame(trendsOut) # convert trendsOut to a dataframe

RR <- exp(trendsOut.df$Cont4.trend) # get the coefficient on the response scale

RR
[1] 1.006696

We can also use the same conversion on the confidence limits of the modelled effect, for the upper:

RR.up <- exp(trendsOut.df$asymp.UCL) # get the upper confidence interval on the response scale
  
RR.up
[1] 1.01052

and lower 95% confidence level:

RR.low <- exp(trendsOut.df$asymp.LCL) # get the lower confidence interval on the response scale

RR.low
[1] 1.002887

This tells you that for a unit change in Cont4, Resp3 increases by 1.007x (95% confidence interval is 1.003 to 1.011) which is an increase of 0.67%.

Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)

This question is only applies to categorical predictors of which their on none in Example 3.

Example 4: Resp4.p ~ Cat5 + Cont6 + Cat5:Cont6 + 1

Recall that Example 4 contains two predictors, one is categorical and one is continuous:

Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1

The best-specified model (now with poisson error distribution assumption) is:

bestMod4.p

Call:  glm(formula = Resp4.p ~ Cat5 + Cont6 + Cat5:Cont6 + 1, family = poisson(link = "log"), 
    data = myDat4.p)

Coefficients:
  (Intercept)        Cat5StB        Cat5StC          Cont6  Cat5StB:Cont6  
     4.374118       0.074861      -0.262595       0.004075      -0.002618  
Cat5StC:Cont6  
     0.003898  

Degrees of Freedom: 49 Total (i.e. Null);  44 Residual
Null Deviance:      183.1 
Residual Deviance: 50.35    AIC: 386.8

that was fit to data in myDat4.p:

str(myDat4.p)
'data.frame':   50 obs. of  3 variables:
 $ Cat5   : Factor w/ 3 levels "StA","StB","StC": 2 1 3 2 2 1 1 3 3 2 ...
 $ Cont6  : num  60.9 55.9 82.8 89.4 53.4 ...
 $ Resp4.p: int  85 112 109 87 99 98 92 101 87 110 ...

and the dredge() table you used to pick your bestMod4.p is

dredgeOut4.p
Global model call: glm(formula = Resp4.p ~ Cat5 + Cont6 + Cat5:Cont6 + 1, family = poisson(link = "log"), 
    data = myDat4.p)
---
Model selection table 
  (Int) Ct5      Cn6 Ct5:Cn6 df   logLik  AICc  delta weight
8 4.374   + 0.004075       +  6 -187.378 388.7   0.00  0.999
4 4.291   + 0.005194          4 -196.412 401.7  13.00  0.001
3 4.261     0.005378          2 -207.753 419.8  31.05  0.000
2 4.669   +                   3 -239.402 485.3  96.62  0.000
1 4.665                       1 -253.742 509.6 120.86  0.000
Models ranked by AICc(x) 

And here’s a visualization of the model effects:

#### i) choosing the values of your predictors at which to make a prediction


# Set up your predictors for the visualized fit
forCat5<-unique(myDat4.p$Cat5) # all levels of your categorical predictor
forCont6<-seq(from = min(myDat4.p$Cont6), to = max(myDat4.p$Cont6), by = 1)# a sequence of numbers in your continuous predictor range
  
# create a data frame with your predictors
forVis<-expand.grid(Cat5 = Cat5, Cont6 = forCont6) # expand.grid() function makes sure you have all combinations of predictors.  

#### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor


# Get your model fit estimates at each value of your predictors
modFit<-predict(bestMod4.p, # the model
                newdata = forVis, # the predictor values
                type = "response", # make the predictions on the response scale
                se.fit = TRUE) # include uncertainty estimate

forVis$Fit<-modFit$fit # add your fit to the data frame
forVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frame
forVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame

#### iii) use the model estimates to plot your model fit


library(ggplot2) # load ggplot2 library

ggplot() + # start ggplot
  
  geom_point(data = myDat4.p, # add observations to your plot
             mapping = aes(x = Cont6, y = Resp4.p, col = Cat5)) + # control position of data points so they are easier to see on the plot
  
  geom_line(data = forVis, # add the modelled fit to your plot
              mapping = aes(x = Cont6, y = Fit, col = Cat5),
              size = 1.2) + # control thickness of line
  
    geom_line(data = forVis, # add uncertainty to your plot (upper line)
              mapping = aes(x = Cont6, y = Upper, col = Cat5),
              size = 0.4, # control thickness of line
              linetype = 2) + # control style of line
  
      geom_line(data = forVis, # add uncertainty to your plot (lower line)
              mapping = aes(x = Cont6, y = Lower, col = Cat5),
              size = 0.4, # control thickness of line
              linetype = 2) + # control style of line
  
  ylab("Resp4, (units)") + # change y-axis label
  
  xlab("Cont6, (units)") + # change x-axis label
  
  labs(fill="Cat5, (units)", col="Cat5, (units)") + # change legend title
  
  theme_bw() # change ggplot theme

What are your modelled effects (with uncertainty)?

The process for estimating effects of categorical vs. continuous predictors is different.

For categorical predictors (Cat5), you can use emmeans() to give you the effects on the response scale directly:

emmOut <- emmeans(object = bestMod4.p, # your model
            specs = ~ Cat5 + Cont6 + Cat5:Cont6 + 1, # your predictors
            type = "response") # report coefficients on the response scale

emmOut
 Cat5 Cont6  rate   SE  df asymp.LCL asymp.UCL
 StA   73.5 107.1 2.32 Inf     102.7       112
 StB   73.5  95.2 2.61 Inf      90.2       100
 StC   73.5 109.7 2.70 Inf     104.5       115

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

When Cat5 is StA, Resp4.p is estimated to be 74 (95% confidence interval: to 103 units).
When Cat5 is StB, Resp4.p is estimated to be 74 (95% confidence interval: to 90 units). When Cat5 is StC, Resp4.p is estimated to be 74 (95% confidence interval: to 105 units).

For continuous predictors (Cont6), you need to use the emtrends() function and then convert the coefficients to the response scale:

trendsOut <- emtrends(bestMod4.p,
                      specs = ~ Cat5 + Cont6 + Cat5:Cont6 + 1, # your predictors
                      var = "Cont6") # your predictor of interest

trendsOut
 Cat5 Cont6 Cont6.trend       SE  df asymp.LCL asymp.UCL
 StA   73.5     0.00408 0.000829 Inf   0.00245   0.00570
 StB   73.5     0.00146 0.001450 Inf  -0.00139   0.00430
 StC   73.5     0.00797 0.000903 Inf   0.00620   0.00974

Confidence level used: 0.95 

Since you have a poisson error distribution assumption, you can convert the estimate made by emtrends() on to the response scale with the exp() function:

trendsOut.df <- data.frame(trendsOut) # convert trendsOut to a dataframe

RR <- exp(trendsOut.df$Cont6.trend) # get the coefficient on the response scale

RR
[1] 1.004084 1.001458 1.008005

You can also use the same conversion on the confidence limits of the modelled effect, for the upper:

RR.up <- exp(trendsOut.df$asymp.UCL) # get the upper confidence interval on the response scale
  
RR.up
[1] 1.005716 1.004309 1.009791

and lower 95% confidence level:

RR.low <- exp(trendsOut.df$asymp.LCL) # get the lower confidence interval on the response scale

RR.low
[1] 1.0024542 0.9986158 1.0062224

Let’s organize these numbers so we can read the effects easier:

RRAll <- data.frame(Cat5 = trendsOut.df$Cat5, # include info on the Cat5 level
                    RR = exp(trendsOut.df$Cont6.trend), # the modelled effect as a rate ratio
                    RR.up = exp(trendsOut.df$asymp.UCL), # upper 95% confidence level of the modelled effect
                    RR.down = exp(trendsOut.df$asymp.LCL)) # lower 95% confidence level of the modelled effect

RRAll
  Cat5       RR    RR.up   RR.down
1  StA 1.004084 1.005716 1.0024542
2  StB 1.001458 1.004309 0.9986158
3  StC 1.008005 1.009791 1.0062224

This tells you that for a unit change in Cont6, Resp4.p

  • increases by 1.004 times (an increase of 0.41%) when Cat5 is StA (95% confidence interval: 1.002 to 1.006).

  • increases by 1.001 times (an increase of 0.15%) when Cat5 is StB (95% confidence interval: 0.999 to 1.004).

  • increases by 1.008 times (an increase of 0.8%) when Cat5 is StC (95% confidence interval: 1.006 to 1.01).

Note that the 95% confidence interval for the effect of Cont6 on Resp4.p (the rate ratio) includes 1. When the rate ratio is 1, there is no effect of the continuous predictor on the response. See Section 5.2.4 for more explanation on this.

You can test for evidence that the effects of Cont6 on Resp4.p are significant with:

test(trendsOut) # get test if coefficient is different than zero.
 Cat5 Cont6 Cont6.trend       SE  df z.ratio p.value
 StA   73.5     0.00408 0.000829 Inf   4.918  <.0001
 StB   73.5     0.00146 0.001450 Inf   1.005  0.3150
 StC   73.5     0.00797 0.000903 Inf   8.828  <.0001

Note that these tests are done on the link scale but can be used for reporting effects on the response scale.

From the results, you can see that:

  • there is evidence that there is an effect of Cont6 on Resp4.p when Cat5 = StA. (i.e. the slope on the link scale is different than zero, and the rate ratio on the response scale is different than 1; P < 0.0001).

  • there is no evidence that there is an effect of Cont6 on Resp4.p when Cat5 = StB. (i.e. the slope on the link scale is not different than zero, and the rate ratio on the response scale is not different than 1; P = 0.32).

  • there is evidence that there is an effect of Cont6 on Resp4.p when Cat5 = StC. (i.e. the slope on the link scale is different than zero, and the rate ratio on the response scale is different than 1; P < 0.0001).

Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)

You can also find out which effects of Cont6 on Resp4.p are different from one another using pairs() on the emtrends() output:

forCompRR <- pairs(trendsOut, 
                   simple = "Cat5")

forCompRR
Cont6 = 73.5:
 contrast  estimate      SE  df z.ratio p.value
 StA - StB  0.00262 0.00167 Inf   1.568  0.2598
 StA - StC -0.00390 0.00123 Inf  -3.180  0.0042
 StB - StC -0.00652 0.00171 Inf  -3.814  0.0004

P value adjustment: tukey method for comparing a family of 3 estimates 

Remember that you need to convert any predictor coefficients to the response scale when reporting.

The table above tells you that:

  • There is no evidence the effect of Cont6 on Resp4.p differs when Cat5 is StA vs. StB (P = 0.26).

  • There is evidence that the effect of Cont6 on Resp4.p is different when Cat5 is StA vs. StC (P = 0.0042). The effect of Cont6 on Resp4.p is exp(-0.0039) = 0.9961076 times as big when Cat5 is StA vs. StC.

  • There is evidence that the effect of Cont6 on Resp4.p is different when Cat5 is StB vs. StC (P = 0.0004). The effect of Cont6 on Resp4.p is exp(-0.00652) = 0.9935012 times as big when Cat5 is StB vs. StC.

Only if all the slopes were similar, you would want to test if the levels of your categorical predictor (Cat5) have significantly different coefficients (intercepts) from each other with pairs() on the output from emmeans() (the emmOut object).


5.2.5 Models with a Gamma error distribution assumption:

Background

For models with a Gamma error distribution assumption, you will typically start with a link function that is the inverse or link = "inverse". You need to take this link function into consideration when you report your modelled effects.

For a model where link = "inverse", the interpretation of the coefficients on the response scale is less clear. In this case, you would use your model to make predictions of your response and report the effects by describing these changes in predictions. We will cover making predictions in the next chapter.

However, you can also use link = "log" for your models with a Gamma error distribution assumption. Using “log” will let you use the reporting methods in the section on “Models with a poisson error distribution assumption”(Section 5.2.4). So, if you have a model with a Gamma error distribution assumption, use link = "log" for modelling, and follow the methods given above (Section 5.2.4).


5.2.6 Models with a binomial error distribution assumption:

Background

For models with a binomial error distribution assumption, you will typically start with a link function that is the logit function or link = "logit". You need to take this link function into consideration when you report your modelled effects.

The logit function is also called the “log-odds” function as it is equal to the logarithm of the odds.

Models with a binomial error distribution assumption model the probability of an event happening. This might be presence (vs. absence), mature (vs. immature), death (vs. alive), etc., but in all cases this is represented as the probability of your response being 1 (vs. 0).

If \(p\) is the probability of an event happening (e.g. response being 1), the probability an event doesn’t happen is \(1-p\) when an event can only have two outcomes. The odds are the probability of the event happening divided by the probability it doesn’t happen, or:

\(\frac{p}{1-p}\)

The logit function (or log-odds function) is then:

\(log_e(\frac{p}{1-p})\)

Jargon alert! Note that a GLM using the logit link function is also called logistic regression.

Alternate link functions for binomially distributed data are probit and cloglog, but the logit link function is by far the one most used as the coefficients that result from the model are the easiest to interpret. But, as always, you should choose the link function that leads to the best-behaved residuals for your model.

Similar to reporting a rate ratio when your link = “log” (Section 5.2.4), you report your coefficients of a model with a binomial error distribution assumption (and link = “logit”) on the response scale by raising e to the power of the coefficients estimated on the link scale.

When you have a continuous predictor, \(e\) is raised to the slope (estimated on the link scale) to give you an odds ratio (OR) on the response scale. 1 - OR gives you the % change in odds for a unit change in your continuous predictor.

When you have a categorical predictor, \(e\) is raised to the intercept of each level of your predictor (estimated on the link scale) to give you the odds associated with that predictor level on the response scale.

Here is the math linking the coefficient to the odds ratio:

Given a binomial model, \[ \begin{align} log_e(\frac{p_i}{1-p_i}) &= \beta_1 \cdot Cov_i + \beta_0 \\ Resp_i &\sim binomial(p_i) \end{align} \] \(p_i\) is the probability of success and \(1-p_i\) is the probability of failure, and the odds are \(\frac{p_i}{1-p_i}\).

Then the odds ratio (OR), or % change in odds due to a change in predictor is:

\[ \begin{align} OR&=\frac{odds\;when\;Cov=a+1}{odds\;when \;Cov=a}\\[2em] OR&=\frac{(\frac{p}{1-p}|Cov=a+1)}{(\frac{p}{1-p}|Cov=a)}\\[2em] log_e(OR)&=(log_e\frac{p}{1-p}|Cov=a+1) - (log_e\frac{p}{1-p}|Cov=a)\\[2em] log_e(OR)&=(\beta_1\cdot (a+1)+\beta_0) - (\beta_1\cdot a+\beta_0)\\[2em] log_e(OR)&=\beta_1\cdot a + \beta_1-\beta_1\cdot a\\[2em] OR &=e^{\beta_1} \end{align} \]

One thing to note:

The emmeans() function will convert the coefficients (intercepts) from the link to the response scale. You can ask for this with the type = "response" argument.

In contrast, the emtrends() function does not convert the coefficients (slopes) to represent effects on the response scale. This is because emtrends() is reporting the slope of a straight line - the trend line on the link scale. But the line isn’t straight on the response scale.

You need to convert from the link to the response by hand when we use emtrends()23. How you make the conversion, and interpret the result, will depend on your error distribution assumption:

Here we will go through examples of reporting for models with a binomial error distribution assumption. If you want more details on why you are using the methods below, check the section on “Models with a normal error distribution assumption” (Section 5.2.3) for context.

Examples

Here I have made new examples where the error distribution assumption is binomial This changes the response variables, indicated by Resp1.b, Resp2.b, etc. to be values of either 0 or 1.

Example 1: Resp1.b ~ Cat1 + 1

Recall that Example 1 contains only one predictor and it is categorical:

Resp1.b ~ Cat1 + 1

The best-specified model (now with a binomial error distribution assumption) is:

bestMod1.b

Call:  glm(formula = Resp1.b ~ Cat1 + 1, family = binomial(link = "logit"), 
    data = myDat1)

Coefficients:
(Intercept)      Cat1StB      Cat1StC  
    -0.6931      -1.8718       1.6094  

Degrees of Freedom: 49 Total (i.e. Null);  47 Residual
Null Deviance:      68.03 
Residual Deviance: 51.43    AIC: 57.43

that was fit to data in myDat1.b:

str(myDat1.b)
'data.frame':   50 obs. of  2 variables:
 $ Cat1   : Factor w/ 3 levels "StA","StB","StC": 3 3 2 3 3 1 1 2 3 3 ...
 $ Resp1.b: int  1 0 0 1 1 0 0 0 1 0 ...

and the dredge() table you used to pick your bestMod1.b is in dredgeOut1.b

dredgeOut1.b
Global model call: glm(formula = Resp1.b ~ Cat1 + 1, family = binomial(link = "logit"), 
    data = myDat1)
---
Model selection table 
  (Intrc) Cat1    R^2 df  logLik AICc delta weight
2 -0.6931    + 0.2825  3 -25.714 57.9  0.00  0.998
1 -0.3228      0.0000  1 -34.015 70.1 12.16  0.002
Models ranked by AICc(x) 

And here’s a visualization of the model effects:

# Set up your predictors for the visualized fit
forCat1<-unique(myDat1.b$Cat1) # every value of your categorical predictor

# create a data frame with your predictors
forVis<-expand.grid(Cat1=forCat1) # expand.grid() function makes sure you have all combinations of predictors

# Get your model fit estimates at each value of your predictors
modFit<-predict(bestMod1.b, # the model
                newdata = forVis, # the predictor values
                type = "response", # make the predictions on the response scale
                se.fit = TRUE) # include uncertainty estimate

forVis$Fit<-modFit$fit # add your fit to the data frame
forVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frame
forVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame


library(ggplot2)

ggplot() + # start ggplot
  geom_point(data = myDat1.b, # add observations to your plot
             mapping = aes(x = Cat1, y = Resp1.b)) + 
  geom_errorbar(data = forVis, # add the uncertainty to your plot
              mapping = aes(x = Cat1, y = Fit, ymin = Lower, ymax = Upper),
              linewidth=1.2) + # control thickness of errorbar line
  geom_point(data = forVis, # add the modelled fit to your plot
              mapping = aes(x = Cat1, y = Fit), 
             shape = 22, # shape of point
             size = 3, # size of point
             fill = "white", # fill colour for plot
             col = 'black') + # outline colour for plot
  ylab("probability Resp1.b = 1")+ # change y-axis label
  xlab("Cat1, (units)")+ # change x-axis label
  theme_bw() # change ggplot theme

What are your modelled effects (with uncertainty)?

For categorical predictors, you can get the coefficients on the response scale directly with the emmeans() function24 when you specify type = “response”. These are the probabilities (not odds or odds ratio) that Resp1.b = 1 at each level of your categorical predictor:

library(emmeans) # load the emmeans package

emmOut <- emmeans(object = bestMod1.b, # your model
            specs = ~ Cat1, # your predictors
            type = "response") # report coefficients on the response scale

emmOut
 Cat1   prob     SE  df asymp.LCL asymp.UCL
 StA  0.3333 0.1220 Inf   0.14596     0.594
 StB  0.0714 0.0688 Inf   0.00996     0.370
 StC  0.7143 0.0986 Inf   0.49239     0.866

Confidence level used: 0.95 
Intervals are back-transformed from the logit scale 

Notice that the column reported is prob: these coefficients are already converted back onto the response scale (as probabilities).

Note also that two types of uncertainty are measured here. SE stands for the standard error around the prediction, and is a measure of uncertainty of the average modelled effect. The lower.CL and upper.CL represent the 95% confidence limits of the prediction - so if you observed a new Resp1.b at a particular Cat1, there is a 95% probability that the probability Resp1.b is 1 would lie within the bounds of the confidence limits.

From this output you see:

When Cat1 is SpA, the probability that Resp1.b is 1 is estimated to be 0.33 (95% confidence interval: 0.15 to 0.59 units). This is equivalent to odds (\(\frac{p}{(1-p)}\)) of 0.5.

When Cat1 is SpB, the probability that Resp1.b is 1 is estimated to be 0.07 (95% confidence interval: 0.01 to 0.37 units). This is equivalent to odds (\(\frac{p}{(1-p)}\)) of 0.08.

When Cat1 is SpC, the probability that Resp1.b is 1 is estimated to be 0.71 (95% confidence interval: 0.49 to 0.87 units). This is equivalent to odds (\(\frac{p}{(1-p)}\)) of 2.5.

Finally, note that you can also get a quick plot of the effects by handing the emmeans() output to plot().

The emmeans package provides an easy way of seeing the effects of Cat1 on Resp1.b the response scale, but you can also get this information from the summary() output:

coef(summary(bestMod1.b))
              Estimate Std. Error   z value   Pr(>|z|)
(Intercept) -0.6931472  0.5477226 -1.265508 0.20568935
Cat1StB     -1.8718022  1.1734223 -1.595165 0.11067536
Cat1StC      1.6094379  0.7302967  2.203814 0.02753745

You can convert the coefficients to odds with the exp() function, for example:

When Cat1 is SpA, the odds that Resp1.b is 1 is estimated to be exp(coef(summary(bestMod1.b))[1]) = 0.5. And you can get the probabilities directly with the plogis() function (plogis(coef(summary(bestMod1.b))[1]) = 0.333).

Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)

Here you can compare the effects of the levels of Cat1 on Resp1.b to see if the effects differ between the levels.

forComp <- pairs(emmOut)

forComp
 contrast  odds.ratio     SE  df null z.ratio p.value
 StA / StB     6.5000 7.6300 Inf    1   1.595  0.2477
 StA / StC     0.2000 0.1460 Inf    1  -2.204  0.0705
 StB / StC     0.0308 0.0352 Inf    1  -3.041  0.0067

P value adjustment: tukey method for comparing a family of 3 estimates 
Tests are performed on the log odds ratio scale 

The output shows the results of the multiple comparison (pairwise) testing.

The values in the p.value column tell you the results of the hypothesis test comparing the coefficients between the two levels of your categorical predictor.

Note the odds.ratio column here. Comparing effects of Cat1 on the probability of Resp1.b is done by comparing odds of Resp1.b being 1 under each level of Cat1.

For example:

  • there is no evidence that the probability of Resp1.b being 1 is different when Cat 1 = SpA vs. when Cat1 = SpB (P = 0.2477).

  • there is little evidence that the probability of Resp1.b being 1 is different when Cat 1 = SpA than when Cat1 = SpC (P = 0.0705). The odds ratio is 0.2.

  • there is strong evidence that the probability of Resp1.b being 1 is different when Cat 1 = SpB than when Cat1 = SpC (P = 0.0067). The odds ratio is 0.0308.

This is a ratio of the odds of that Resp1.b is 1 when Cat 1 = SpA vs. the odds of that Resp1.b is 1 when Cat 1 = SpC:

For example, the odds that Resp1.b is 1 when Cat 1 = SpA are:

pCat1.SpA <- summary(emmOut)[1,2]

oddsCat1.SpA <- pCat1.SpA /(1-pCat1.SpA)

oddsCat1.SpA
[1] 0.5

and the odds that Resp1.b is 1 when Cat 1 = SpC are:

pCat1.SpC <- summary(emmOut)[3,2]

oddsCat1.SpC <- pCat1.SpC /(1-pCat1.SpC)

oddsCat1.SpC
[1] 2.5

and the odds ratio StA/StB is

oddsCat1.SpA/oddsCat1.SpC
[1] 0.2

Note that you can also get the results from the pairwise testing visually by handing the output of pairs() to plot().

Example 2: Resp2.b ~ Cat2 + Cat3 + Cat2:Cat3 + 1

Recall that Example 2 contains two predictors and both are categorical:

Resp2.b ~ Cat2 + Cat3 + Cat2:Cat3 + 1

The best-specified model (now with binomial error distribution assumption) is:

bestMod2.b

Call:  glm(formula = Resp2.b ~ Cat2 + Cat3 + Cat2:Cat3 + 1, family = binomial(link = "logit"), 
    data = myDat2.b)

Coefficients:
      (Intercept)            Cat2StB            Cat2StC          Cat3Treat  
         -0.15415            0.23639           -0.08352           -0.10368  
Cat2StB:Cat3Treat  Cat2StC:Cat3Treat  
         -0.29222            1.99958  

Degrees of Freedom: 499 Total (i.e. Null);  494 Residual
Null Deviance:      692.9 
Residual Deviance: 649.7    AIC: 661.7

that was fit to data in myDat2.b:

str(myDat2.b)
'data.frame':   500 obs. of  3 variables:
 $ Cat2   : Factor w/ 3 levels "StA","StB","StC": 2 2 3 2 2 3 2 1 1 1 ...
 $ Cat3   : Factor w/ 2 levels "Control","Treat": 1 1 1 2 1 2 1 2 2 1 ...
 $ Resp2.b: int  0 0 0 1 1 1 0 0 0 1 ...

and the dredge() table you used to pick your bestMod2.b is in

dredgeOut2.b
Global model call: glm(formula = Resp2.b ~ Cat2 + Cat3 + Cat2:Cat3 + 1, family = binomial(link = "logit"), 
    data = myDat2.b)
---
Model selection table 
     (Int) Ct2 Ct3 Ct2:Ct3      R^2 df   logLik  AICc delta weight
8 -0.15420   +   +       + 0.082720  6 -324.843 661.9  0.00      1
4 -0.38060   +   +         0.031870  4 -338.331 684.7 22.89      0
2 -0.20190   +             0.023290  3 -340.538 687.1 25.27      0
3 -0.11690       +         0.007163  2 -344.632 693.3 31.43      0
1  0.04801                 0.000000  1 -346.430 694.9 33.01      0
Models ranked by AICc(x) 

And here’s a visualization of the model effects:

#### i) choosing the values of your predictors at which to make a prediction

# Set up your predictors for the visualized fit
forCat2<-unique(myDat2.b$Cat2) # every level of your categorical predictor
forCat3<-unique(myDat2.b$Cat3) # every level of your categorical predictor
  
# create a data frame with your predictors
forVis<-expand.grid(Cat2 = forCat2, Cat3 = forCat3) # expand.grid() function makes sure you have all combinations of predictors

#### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor


# Get your model fit estimates at each value of your predictors
modFit<-predict(bestMod2.b, # the model
                newdata = forVis, # the predictor values
                type = "response", # make the predictions on the response scale
                se.fit = TRUE) # include uncertainty estimate

forVis$Fit<-modFit$fit # add your fit to the data frame
forVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frame
forVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame

#### iii) use the model estimates to plot your model fit

library(ggplot2) # load ggplot2 library

ggplot() + # start ggplot
  geom_point(data = myDat2.b, # add observations to your plot
             mapping = aes(x = Cat2, y = Resp2.b, col = Cat3), 
             position=position_jitterdodge(jitter.width=0.75, dodge.width=0.75)) + # control position of data points so they are easier to see on the plot
  geom_errorbar(data = forVis, # add the uncertainty to your plot
              mapping = aes(x = Cat2, y = Fit, ymin = Lower, ymax = Upper, col = Cat3),
              position=position_dodge(width=0.75), # control position of data points so they are easier to see on the plot
              size=1.2) + # control thickness of errorbar line
  geom_point(data = forVis, # add the modelled fit to your plot
              mapping = aes(x = Cat2, y = Fit, fill = Cat3), 
             shape = 22, # shape of point
             size=3, # size of point
             col = 'black', # outline colour for point
             position=position_dodge(width=0.75)) + # control position of data points so they are easier to see on the plot
  ylab("probability Resp2.b = 1")+ # change y-axis label
  xlab("Cat2, (units)")+ # change x-axis label
  labs(fill="Cat3, (units)", col="Cat3, (units)") + # change legend title
  theme_bw() # change ggplot theme

What are your modelled effects (with uncertainty)?

As above, we can use the emmeans package to see the effects of the predictors on the scale of the response:

emmOut <- emmeans(object = bestMod2.b, # your model
            specs = ~ Cat2 + Cat3 + Cat2:Cat3, # your predictors
            type = "response") # report coefficients on the response scale

emmOut
 Cat2 Cat3     prob     SE  df asymp.LCL asymp.UCL
 StA  Control 0.462 0.0523 Inf     0.362     0.564
 StB  Control 0.521 0.0585 Inf     0.407     0.632
 StC  Control 0.441 0.0515 Inf     0.344     0.543
 StA  Treat   0.436 0.0561 Inf     0.331     0.547
 StB  Treat   0.422 0.0521 Inf     0.325     0.526
 StC  Treat   0.840 0.0423 Inf     0.739     0.907

Confidence level used: 0.95 
Intervals are back-transformed from the logit scale 

So now you have a modelled value of your response for each level of your categorical predictors. For example:

When Cat2 is StA and Cat3 is Treat, the probability that Resp1.b is 1 is estimated to be 0.44 (95% confidence interval: 0.33 to 0.55 units).

When Cat2 is StB and Cat3 is Control, the probability that Resp1.b is 1 is estimated to be 0.52 (95% confidence interval: 0.41 to 0.63 units).

Finally, note that you can also get a quick plot of the effects by handing the emmeans() output to plot().

Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)

As with Example 1, you can also find out which combinations of predictor levels are leading to statistically different model predictions in Resp2.b:

forComp <- pairs(emmOut, # emmeans output
                 simple = "Cat2") # contrasts within this categorical predictor

forComp
Cat3 = Control:
 contrast  odds.ratio     SE  df null z.ratio p.value
 StA / StB      0.789 0.2490 Inf    1  -0.751  0.7331
 StA / StC      1.087 0.3220 Inf    1   0.282  0.9572
 StB / StC      1.377 0.4320 Inf    1   1.019  0.5646

Cat3 = Treat:
 contrast  odds.ratio     SE  df null z.ratio p.value
 StA / StB      1.057 0.3300 Inf    1   0.179  0.9826
 StA / StC      0.147 0.0573 Inf    1  -4.925  <.0001
 StB / StC      0.139 0.0530 Inf    1  -5.183  <.0001

P value adjustment: tukey method for comparing a family of 3 estimates 
Tests are performed on the log odds ratio scale 

Based on a threshold P-value of 0.05, you can see that the effects at some combinations of your predictors are not statistically different from each other. For example, there is no evidence that the probability of Resp2.b being 1 when Cat3 = Control and Cat2 = StA (0.46)25 than when Cat3 = Control and Cat2 = StB (0.52) are different from each other (P = 0.73).

On the other hand, some combinations of your predictors are statistically different from each other. For example, there is strong evidence that the probability of Resp2.b being 1 when Cat3 = Treat and Cat2 = StA (0.46) is less than that when Cat3 = Treat and Cat2 = StC (0.84) are different from each other (P < 0.0001).

Note that you can also get the results from the pairwise testing visually by handing the output of pairs() to plot().

Example 3: Resp3.b ~ Cont4 + 1

Recall that Example 3 contains one predictor and the predictor is continuous:

Resp3 ~ Cont4 + 1

The best-specified model (now with a binomial error distribution assumption) is:

bestMod3.b

Call:  glm(formula = Resp3.b ~ Cont4 + 1, family = binomial(link = "logit"), 
    data = myDat3.b)

Coefficients:
(Intercept)        Cont4  
    -8.1519       0.1086  

Degrees of Freedom: 299 Total (i.e. Null);  298 Residual
Null Deviance:      415.8 
Residual Deviance: 199.1    AIC: 203.1

that was fit to data in myDat3.b:

str(myDat3.b)
'data.frame':   300 obs. of  2 variables:
 $ Cont4  : num  57.7 53.2 79.7 35.1 72.2 ...
 $ Resp3.b: int  1 0 1 0 0 1 1 1 1 0 ...

and the dredge() table you used to pick your bestMod3.b is

dredgeOut3.b
Global model call: glm(formula = Resp3.b ~ Cont4 + 1, family = binomial(link = "logit"), 
    data = myDat3.b)
---
Model selection table 
   (Intrc)  Cont4    R^2 df   logLik  AICc  delta weight
2 -8.15200 0.1086 0.5144  2  -99.533 203.1   0.00      1
1  0.04001        0.0000  1 -207.884 417.8 214.68      0
Models ranked by AICc(x) 

And here’s a visualization of the model effects:

#### i) choosing the values of your predictors at which to make a prediction


# Set up your predictors for the visualized fit
forCont4<-seq(from = min(myDat3.b$Cont4), to = max(myDat3.b$Cont4), by = 1)# a sequence of numbers in your continuous predictor range
  
# create a data frame with your predictors
forVis<-expand.grid(Cont4 = forCont4) # expand.grid() function makes sure you have all combinations of predictors.  

#### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor


# Get your model fit estimates at each value of your predictors
modFit<-predict(bestMod3.b, # the model
                newdata = forVis, # the predictor values
                type = "response", # make the predictions on the response scale
                se.fit = TRUE) # include uncertainty estimate

forVis$Fit<-modFit$fit # add your fit to the data frame
forVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frame
forVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame

#### iii) use the model estimates to plot your model fit


library(ggplot2) # load ggplot2 library

ggplot() + # start ggplot
  
  geom_point(data = myDat3.b, # add observations to your plot
             mapping = aes(x = Cont4, y = Resp3.b)) + # control position of data points so they are easier to see on the plot
  
  geom_line(data = forVis, # add the modelled fit to your plot
              mapping = aes(x = Cont4, y = Fit),
              size = 1.2) + # control thickness of line
  
    geom_line(data = forVis, # add uncertainty to your plot (upper line)
              mapping = aes(x = Cont4, y = Upper),
              size = 0.8, # control thickness of line
              linetype = 2) + # control style of line
  
      geom_line(data = forVis, # add uncertainty to your plot (lower line)
              mapping = aes(x = Cont4, y = Lower),
              size = 0.8, # control thickness of line
              linetype = 2) + # control style of line
  
  ylab("probability Resp2.b = 1")+ # change y-axis label
  
  xlab("Cont4, (units)") + # change x-axis label
  
  theme_bw() # change ggplot theme

The model fit line may seem a little strange as the data are binary (can only be 0 or 1) but the model fit line goes through the space in between 0 and 1. This model fit line tells you how the probability of Resp3.b being 1 changes as Cont4 changes.

What are your modelled effects (with uncertainty)?

trendsOut <- emtrends(bestMod3.b,
                      specs = ~ Cont4 + 1, # your predictors
                      var = "Cont4") # your predictor of interest

trendsOut
 Cont4 Cont4.trend     SE  df asymp.LCL asymp.UCL
  75.8       0.109 0.0116 Inf    0.0859     0.131

Confidence level used: 0.95 

Note that this is the same estimate as you find in the summary() output of your model object:

coef(summary(bestMod3.b))
              Estimate Std. Error   z value     Pr(>|z|)
(Intercept) -8.1518996 0.89160247 -9.142976 6.075711e-20
Cont4        0.1085946 0.01160302  9.359159 8.037377e-21

This is because emtrends() returns coefficients on the link scale.

To report the coefficients on a response scale, you need to consider your link function. Since you have a binomial error distribution assumption, you can convert the estimate made by emtrends() to the response scale with:

trendsOut.df <- data.frame(trendsOut) # convert trendsOut to a dataframe

OR <- exp(trendsOut.df$Cont4.trend) # get the coefficient on the response scale

OR
[1] 1.11471

You can also use the same conversion on the confidence limits of the modelled effect, for the upper:

OR.up <- exp(trendsOut.df$asymp.UCL) # get the upper confidence interval on the response scale
  
OR.up
[1] 1.140351

and lower 95% confidence level:

OR.low <- exp(trendsOut.df$asymp.LCL) # get the lower confidence interval on the response scale

OR.low
[1] 1.089646

This coefficient - called the odds ratio (Section 5.2.6) - tells you that for a unit change in Cont4, the odds that Resp3.b is 1 increases by 1.115 times (95% confidence interval is 1.09 to 1.14) which is an increase in the odds of Resp3.b being 1 of 11.5%.

You can get a better sense of the odds ratio (vs. odds vs. probabilities) by estimating it directly from the modelled probability that Resp3.b is 1:

Let’s find the probability of Resp3.b = 1 when Cont4 is 60:

pCont4.60 <- predict(bestMod3.b, # model
                 newdata = data.frame(Cont4 = 60), # predictor values for the prediction
                 type = "response", # give prediction on the response scale
                 se.fit = TRUE) # include uncertainty

pCont4.60 <- pCont4.60$fit

pCont4.60
        1 
0.1629792 

So there is a 16.3% probability of Resp3.b being 1 when Cont4 is 60.

The odds associated with Resp3.b being 1 when Cont4 is 60 is given by the probability of Cont4 is divided by the probability of absence:

oddsCont4.60 <- pCont4.60/(1-pCont4.60)

oddsCont4.60
        1 
0.1947134 

Now let’s find the probability of Resp3.b being 1 when Cont4 is 61:

pCont4.61 <- predict(bestMod3.b, # model
                 newdata = data.frame(Cont4 = 61), # predictor values for the prediction
                 type = "response", # give prediction on the response scale
                 se.fit = TRUE) # include uncertainty

pCont4.61 <- pCont4.61$fit

pCont4.61
        1 
0.1783404 

So there is a 17.8% probability of Resp3.b being 1 when Cont4 is 61.

The odds associated with Resp3.b being 1 when Cont4 is 61 is given by the probability of presence divided by the probability of absence:

oddsCont4.61 <- pCont4.61/(1-pCont4.61)

oddsCont4.61
       1 
0.217049 

Now the odds ratio is found by comparing the change in odds when your predictor (Cont4) changes one unit

oddsCont4.61/oddsCont4.60
      1 
1.11471 

Note that this is identical to the calculation of the odds ratio from the coefficient above:

OR
[1] 1.11471

and you interpret this as a OR - 1 or 11.5% change in the odds of Resp3.b being 1 with a unit change of Cont4.

Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)

This question only applies to categorical predictors of which their are none in Example 3.

Example 4: Resp4.b ~ Cat5 + Cont6 + Cat5:Cont6 + 1

Recall that Example 4 contains two predictors, one is categorical and one is continuous:

Resp4.b ~ Cat5 + Cont6 + Cat5:Cont6 + 1

The best-specified model (now with a binomial error distribution assumption) is:

bestMod4.b

Call:  glm(formula = Resp4.b ~ Cat5 + Cont6 + Cat5:Cont6 + 1, family = binomial(link = "logit"), 
    data = myDat4.b)

Coefficients:
  (Intercept)        Cat5StB        Cat5StC          Cont6  Cat5StB:Cont6  
    -28.55811       19.23667      -10.40094        0.15899       -0.09544  
Cat5StC:Cont6  
      0.03623  

Degrees of Freedom: 499 Total (i.e. Null);  494 Residual
Null Deviance:      639.7 
Residual Deviance: 169  AIC: 181

that was fit to data in myDat4.b:

str(myDat4.b)
'data.frame':   500 obs. of  3 variables:
 $ Cat5   : Factor w/ 3 levels "StA","StB","StC": 2 1 1 1 2 1 2 1 3 2 ...
 $ Cont6  : num  206 191 312 223 210 ...
 $ Resp4.b: int  1 1 1 1 1 1 1 0 1 1 ...

and the dredge() table you used to pick your bestMod4.b is

dredgeOut4.b
Global model call: glm(formula = Resp4.b ~ Cat5 + Cont6 + Cat5:Cont6 + 1, family = binomial(link = "logit"), 
    data = myDat4.b)
---
Model selection table 
     (Int) Ct5     Cn6 Ct5:Cn6     R^2 df   logLik  AICc  delta weight
8 -28.5600   + 0.15900       + 0.60990  6  -84.513 181.2   0.00      1
4 -19.4300   + 0.10820         0.59390  4  -94.548 197.2  15.98      0
3 -12.4100     0.07056         0.49520  2 -148.970 302.0 120.77      0
2   0.5276   +                 0.07603  3 -300.081 606.2 425.01      0
1   0.6722                     0.00000  1 -319.850 641.7 460.51      0
Models ranked by AICc(x) 

And here’s a visualization of the model effects:

#### i) choosing the values of your predictors at which to make a prediction


# Set up your predictors for the visualized fit
forCat5<-unique(myDat4.b$Cat5) # all levels of your categorical predictor
forCont6<-seq(from = min(myDat4.b$Cont6), to = max(myDat4.b$Cont6), by = 1)# a sequence of numbers in your continuous predictor range
  
# create a data frame with your predictors
forVis<-expand.grid(Cat5 = forCat5, Cont6 = forCont6) # expand.grid() function makes sure you have all combinations of predictors.  

#### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor


# Get your model fit estimates at each value of your predictors
modFit<-predict(bestMod4.b, # the model
                newdata = forVis, # the predictor values
                type = "response", # make the predictions on the response scale
                se.fit = TRUE) # include uncertainty estimate

forVis$Fit<-modFit$fit # add your fit to the data frame
forVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frame
forVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame

#### iii) use the model estimates to plot your model fit


library(ggplot2) # load ggplot2 library

ggplot() + # start ggplot
  
  geom_point(data = myDat4.b, # add observations to your plot
             mapping = aes(x = Cont6, y = Resp4.b, col = Cat5)) + # control position of data points so they are easier to see on the plot
  
  geom_line(data = forVis, # add the modelled fit to your plot
              mapping = aes(x = Cont6, y = Fit, col = Cat5),
              size = 1.2) + # control thickness of line
  
    geom_line(data = forVis, # add uncertainty to your plot (upper line)
              mapping = aes(x = Cont6, y = Upper, col = Cat5),
              size = 0.4, # control thickness of line
              linetype = 2) + # control style of line
  
      geom_line(data = forVis, # add uncertainty to your plot (lower line)
              mapping = aes(x = Cont6, y = Lower, col = Cat5),
              size = 0.4, # control thickness of line
              linetype = 2) + # control style of line
  
  ylab("probability Resp2.b = 1")+ # change y-axis label
  
  xlab("Cont6, (units)") + # change x-axis label
  
  labs(fill="Cat5, (units)", col="Cat5, (units)") + # change legend title
  
  theme_bw() # change ggplot theme

What are your modelled effects (with uncertainty)?

The process for estimating effects of categorical vs. continuous predictors is different.

For categorical predictors (Cat5), you can use emmeans() to give you the effects on the response scale directly:

emmOut <- emmeans(object = bestMod4.b, # your model
            specs = ~ Cat5 + Cont6 + Cat5:Cont6 + 1, # your predictors
            type = "response") # report coefficients on the response scale

emmOut
 Cat5 Cont6  prob     SE  df asymp.LCL asymp.UCL
 StA    198 0.951 0.0346 Inf     0.819     0.988
 StB    198 0.964 0.0178 Inf     0.907     0.986
 StC    198 0.438 0.0941 Inf     0.269     0.622

Confidence level used: 0.95 
Intervals are back-transformed from the logit scale 

This output tells you that:

  • When Cat5 is StA and Cont6 is 227 units26, there is a 95% probability that Resp4.b will be 1 (95% confidence interval: 0.82 to 0.99%)

  • When Cat5 is StB and Cont6 is 227 units27, there is a 96% probability that Resp4.b will be 1 (95% confidence interval: 0.91 to 0.99 units).

  • When Cat5 is StB and Cont6 is 227 units28, there is a 44% probability that Resp4.b will be 1 (95% confidence interval: 0.27 to 0.62 units).

Note that emmeans() sets your continuous predictor (Cont6) to the mean value of Cont6 (227 units). You can also control this with the at = argument in the emmeans() function:

emmOut <- emmeans(object = bestMod4.b, # your model
            specs = ~ Cat5 + Cont6 + Cat5:Cont6 + 1, # your predictors
            type = "response", # report coefficients on the response scale
            at = list(Cont6 = 150)) # control the value of your continuous predictor at which to make the coefficient estimates


emmOut
 Cat5 Cont6     prob       SE  df asymp.LCL asymp.UCL
 StA    150 8.93e-03 0.009230 Inf  1.17e-03   0.06500
 StB    150 5.53e-01 0.078500 Inf  3.99e-01   0.69695
 StC    150 6.27e-05 0.000142 Inf  7.00e-07   0.00524

Confidence level used: 0.95 
Intervals are back-transformed from the logit scale 

For continuous predictors (Cont6), you need to use the emtrends() function and then convert the coefficients to the response scale:

trendsOut <- emtrends(bestMod4.b,
                      specs = ~ Cat5 + Cont6 + Cat5:Cont6 + 1 , # your predictors
                      var = "Cont6") # your predictor of interest

trendsOut
 Cat5 Cont6 Cont6.trend     SE  df asymp.LCL asymp.UCL
 StA    198      0.1590 0.0330 Inf    0.0942    0.2237
 StB    198      0.0635 0.0114 Inf    0.0413    0.0858
 StC    198      0.1952 0.0444 Inf    0.1083    0.2821

Confidence level used: 0.95 

Since you have a binomial error distribution assumption, you can convert the estimate made by emtrends() on to the response scale with the exp() function:

trendsOut.df <- data.frame(trendsOut) # convert trendsOut to a dataframe

OR <- exp(trendsOut.df$Cont6.trend) # get the coefficient on the response scale

OR
[1] 1.172328 1.065611 1.215576

You can also use the same conversion on the confidence limits of the modelled effect, for the upper:

OR.up <- exp(trendsOut.df$asymp.UCL) # get the upper confidence interval on the response scale
  
OR.up
[1] 1.250737 1.089631 1.325976

and lower 95% confidence level:

OR.low <- exp(trendsOut.df$asymp.LCL) # get the lower confidence interval on the response scale

OR.low
[1] 1.098834 1.042120 1.114368

Let’s organize these numbers so we can read the effects easier:

ORAll <- data.frame(Cat5 = trendsOut.df$Cat5, # include info on the Cat5 level
                    OR = exp(trendsOut.df$Cont6.trend), # the modelled effect as a rate ratio
                    OR.up = exp(trendsOut.df$asymp.UCL), # upper 95% confidence level of the modelled effect
                    OR.down = exp(trendsOut.df$asymp.LCL)) # lower 95% confidence level of the modelled effect

ORAll
  Cat5       OR    OR.up  OR.down
1  StA 1.172328 1.250737 1.098834
2  StB 1.065611 1.089631 1.042120
3  StC 1.215576 1.325976 1.114368

This tells you that for a unit change in Cont6, the odds that Resp4.b is 1:

  • increases by 17.2% with a unit change in Cont6 when Cat5 is StA (95% confidence interval: 9.9 to 25.1).

  • increases by 6.6% with a unit change in Cont6 when Cat5 is StB (95% confidence interval: 4.2 to 9).

  • increases by 21.6% with a unit change in Cont6 when Cat5 is StC (95% confidence interval: 11.4 to 32.6).

You can test for evidence that the effects of Cont6 on Resp4.b are significant (different than zero) with:

test(trendsOut) # get test if coefficient is different than zero.
 Cat5 Cont6 Cont6.trend     SE  df z.ratio p.value
 StA    198      0.1590 0.0330 Inf   4.813  <.0001
 StB    198      0.0635 0.0114 Inf   5.588  <.0001
 StC    198      0.1952 0.0444 Inf   4.401  <.0001

Note that these tests are done on the link scale but can be used for reporting effects on the response scale.

From the results, you can see that the effects of Cont6 on Resp4.b at all levels of Cat5 are significant (P < 0.0001).

Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)

You can also find out which effects of Cont6 on Resp4.b are different from one another using pairs() on the emtrends() output:

forCompOR <- pairs(trendsOut, 
                   simple = "Cat5")

forCompOR
Cont6 = 198:
 contrast  estimate     SE  df z.ratio p.value
 StA - StB   0.0954 0.0349 Inf   2.732  0.0173
 StA - StC  -0.0362 0.0553 Inf  -0.655  0.7895
 StB - StC  -0.1317 0.0458 Inf  -2.876  0.0112

P value adjustment: tukey method for comparing a family of 3 estimates 

The table above tells you that:

  • There is evidence that the effect of Cont6 on Resp4.b is different when Cat5 is StA vs. StB (P = 0.017). The effect of Cont6 on Resp4.b is exp(0.0954) = 1.1 times as big when Cat5 is StA vs. StB (i.e. the effect of Cont6 on Resp4.b is 10% more when Cat5 is StA vs. StB)

  • There is no evidence the effect of Cont6 on Resp4.b differs when Cat5 is StA vs. StC (P = 0.789).

  • There is evidence that the effect of Cont6 on Resp4.b is different when Cat5 is StB vs. StC (P = 0.011). The effect of Cont6 on Resp4.b is exp(-0.1317) = 0.88 times as big when Cat5 is StB vs. StC (i.e. the effect of Cont6 on Resp4.b is 12% less when Cat5 is StB vs. StC)

Only if all the slopes29 were similar, you would want to test if the levels of your categorical predictor (Cat5) have significantly different coefficients (intercepts) from each other with pairs() on the output from emmeans() (the emmOut object).

Up next

In the next chapter, we will discuss how you can use your model to make predictions.

6 Communicating your hypothesis testing results

Here we’ll walk through your statistical modelling framework and discuss how to communicate motivations, methods and results of your data analysis and hypothesis testing for a biological report/article (& the exam!).

Here again is our statistical modelling framework:

Note that I’ve merged the Visualizing section with the Reporting section to correspond with how we will write-up our results.

6.1 A familiar example

I’ll use the following example case as I discuss how to communicate your results:

You are interested in moose (or elk, Alces alces) and how their growth rate changes from population to population. You find information on mean annual weight change by moose in 12 locations in Sweden. You wonder if latitude can explain some of the variability in growth rates you are observing. You also wonder if minimum winter temperature might explain some variability in growth rate. Finally, you wonder if changes in growth rate with latitude and/or winter temperature depend on each other (i.e, latitudinal dependent effects on growth rate depend on what the minimum winter temperature is) and/or sex (male and female) - i.e. the environmentally dependent change in growth rate depends on what sex the animal is.


And here’s a screenshot of the readme file for the moose data:

6.2 Response(s)

The task:

Your response variable30 is the observed variability you are trying to explain. Your response forms your research question - i.e. “why is my response varying?”. In this section, you need to identify the variability you’re trying to explain and why (your motivation).

The code:

Code supporting your tasks in this section allow you to import your data and help you describe your response variable, e.g.:

rm(list=ls()) # make sure your workspace is clear

#Import the data
myDat <- read.table("ExDat.csv", # file name
                    sep = ',',  # lines in the file are separated by a comma - this information is in the README file
                    dec = '.', # decimals are denoted by '.' - this information is in the README file
                    stringsAsFactors = TRUE, # make sure strings (characters) are imported as factors
                    header = TRUE) # there is a column header in the file

str(myDat) # check the structure of the data file
'data.frame':   24 obs. of  6 variables:
 $ Study  : int  1 2 3 4 5 6 7 8 9 12 ...
 $ Sex    : Factor w/ 2 levels "Females","Males": 1 1 1 1 1 1 1 1 1 1 ...
 $ delMass: num  40.2 38.7 41.2 39.4 36.2 40 44.1 42.2 46.2 43.6 ...
 $ Lat    : num  58 57.7 58 57.9 59.8 61.6 62 62.7 64 65.5 ...
 $ WTemp  : num  -2.4 -1.7 -1.7 -1.7 -4.2 -5.6 -8.1 -7.7 -7.6 -9.3 ...
 $ Source : Factor w/ 3 levels "Gothenburg","Lund",..: 2 3 3 2 3 2 1 1 3 3 ...
summary(myDat) # get summary information of each column
     Study            Sex        delMass           Lat            WTemp        
 Min.   : 1.00   Females:12   Min.   :36.20   Min.   :57.70   Min.   :-11.900  
 1st Qu.: 3.75   Males  :12   1st Qu.:41.95   1st Qu.:58.00   1st Qu.: -8.400  
 Median : 6.50                Median :47.35   Median :61.80   Median : -6.600  
 Mean   : 7.00                Mean   :47.77   Mean   :61.60   Mean   : -6.100  
 3rd Qu.: 9.75                3rd Qu.:54.23   3rd Qu.:64.38   3rd Qu.: -2.225  
 Max.   :14.00                Max.   :60.60   Max.   :66.00   Max.   : -1.700  
        Source  
 Gothenburg: 6  
 Lund      : 8  
 Stockholm :10  
                
                
                

The write-up:

Use the code above, the data readme file, as well as supporting literature to communicate about your response variable and motivation. Include:

  • what your response variable is

  • why you want to explain variability in your response

  • your research question

  • how your response variable was measured and a description of how your response varies. For example, give the range in your response if possible. If your response is a category, give the levels of the category (e.g. alive vs. dead). Always give units, if applicable.

For example:

Here I want to explain variability in growth rate changes in moose (Alces alces). Moose size will be related to survival rates, and likely also their ability to reproduce (e.g. Sand et al. 1995). Being able to explain what controls growth rate in moose may help us explain why some populations are more productive than other populations, and it will help us predict productivity for populations at other places or times.

My research question is “why does growth rate vary?”. My observations of growth rate were measured from adult moose across 12 locations in Sweden. The growth rate varied from 36.2 to 60.6 kg per year.

6.3 Covariate(s)

The task:

Here you will identify the covariates in your model - i.e. what variables might explain variability in your response. To write-up this section, start with the biological mechanism (or processes) that make you think the covariate could be influencing your response. Then describe how that covariate was measured. At this point you might identify mechanisms that you are not able to measure, e.g. because they were beyond the scope of your study. Make note of these. You will have a chance to talk about them at the end of this exercise (in Reporting).

The code:

Code supporting your tasks in this section will help you describe your covariates, e.g.:

# Getting information on my covariates
range(myDat$Lat) # ranges for continuous covariates
[1] 57.7 66.0
range(myDat$WTemp) # ranges for continuous covariates
[1] -11.9  -1.7
unique(myDat$Sex) # factor levels for categorical covariates
[1] Females Males  
Levels: Females Males

The write-up:

Use the code above, the data readme file, as well as supporting literature to communicate about your covariate variables. Include:

  • what your covariate variables are

  • why you think each covariate might explain variability in your response (mechanisms)

  • how each covariate variable was measured and a description of the covariate. This means at least giving the range of each covariate (for continuous covariates) or the levels of any categorical covariate. Remember to include units where applicable.

For example:

Here I consider the effects of latitude (˚N), minimum winter temperature (˚C) and sex (male vs. female) on annual growth in adult moose.

Latitude may impact the growth of the moose as the environment changes with latitude (e.g. temperature, food, light level, growing season). Differences in food availability may lead to growth differences [citation]31. In this study, latitude is measured in decimal degrees N and ranges from 57.7 to 66.0˚N.

Winter harshness may impact annual moose growth by making energetic losses larger in the winter due to increased costs of temperature regulation [citation]. Here, I measure winter harshness as minimum winter temperature (˚C) which ranged from -11.9 to -1.7˚C in my study.

Growth differences might also occur due to sex of the animal as male and female moose have different life history strategies (e.g. reproductive investment)[citation]. Here, I include effects of sex with observations of male vs. female for each growth rate measurement.

6.4 Hypothesis

The task:

In this section you will describe your research hypothesis including main effects and any interactions. Remember to include the formula notation we’ve been practicing in class (i.e. Response ~ Covariate…)

The code:

No code needed here!

The write-up:

Use your answers to the previous sections to state your research hypothesis. Remember to consider interactions as well as the main effects. Include:

  • a description in words

  • a description in the formula notation

  • a definition of each term

For example,

I will test the research hypothesis that varibility in adult moose growth rate (delMass, kg per year) is explained by latitude (Lat, ˚N), minimum winter temperature (WTemp, ˚C) and sex (Sex, male or female). I will test this by modelling:

delMass ~ Lat + WTemp + Sex + Lat:Sex + WTemp:Sex + Lat:WTemp + Lat:WTemp:Sex

which includes main effects and all possible interaction effects of my covariates (where “:” indicates an interaction between covariates).

6.5 Starting model

The task:

This step involves describing how you built the statistical model to test your hypothesis. Here, you will:

  • Choose your error distribution assumption

  • Choose your shape assumption

  • Choose your starting model

  • Fit your starting model

The code:

Code supporting your tasks in this section will help you support your choices of your model assumptions, e.g.:

# Choosing your error distribution assumption:
## Understanding the nature of your response variable
summary(myDat$delMass) # a summary description of delMass
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  36.20   41.95   47.35   47.77   54.23   60.60 
# Choosing your shape assumption:
## This can include plotting your response variable vs. each continuous covariate
### For Lat:
library(ggplot2) # load ggplot2 library
ggplot()+ # start ggplot
  geom_point(data = myDat,
             mapping = aes(x = Lat, y = delMass))+ # add observations
  xlab("Latitude, ˚N")+ # change x-axis label
  ylab("Annual growth rate, kg per year")+ # change y-axis label
  labs(caption = "Figure 1: Observed annual growth rate (kg per year) vs. latitude (˚N)")+ # add figure caption
  theme_bw()+ # change theme of plot
  theme(plot.caption = element_text(hjust=0)) # move figure legend (caption) to left alignment. Use hjust = 0.5 to align in the center.

#### Relationship between Lat and delMass looks linear


### For WTemp:
ggplot()+ # start ggplot
  geom_point(data = myDat,
             mapping = aes(x = WTemp, y = delMass))+ # add observations
  xlab("Minimum winter temperature, ˚C")+ # change x-axis label
  ylab("Annual growth rate, kg per year")+ # change y-axis label
  labs(caption = "Figure 2: Observed annual growth rate (kg per year) vs. minimum winter temperature (˚C)")+ # add figure caption
  theme_bw()+ # change theme of plot
  theme(plot.caption = element_text(hjust=0)) # move figure legend (caption) to left alignment. Use hjust = 0.5 to align in the center.

#### Relationship between WTemp and delMass looks linear


# Fitting your starting model
startMod<-glm(formula = delMass ~ Lat + WTemp + Sex + Lat:Sex + WTemp:Sex + Lat:WTemp + Lat:WTemp:Sex, # hypothesis
              data = myDat, # data
              family = gaussian(link = "identity")) # error distribution assumption

startMod

Call:  glm(formula = delMass ~ Lat + WTemp + Sex + Lat:Sex + WTemp:Sex + 
    Lat:WTemp + Lat:WTemp:Sex, family = gaussian(link = "identity"), 
    data = myDat)

Coefficients:
       (Intercept)                 Lat               WTemp            SexMales  
          13.86368             0.43887             1.37886           -54.84374  
      Lat:SexMales      WTemp:SexMales           Lat:WTemp  Lat:WTemp:SexMales  
           1.09654            -5.46475            -0.02401             0.08799  

Degrees of Freedom: 23 Total (i.e. Null);  16 Residual
Null Deviance:      1200 
Residual Deviance: 75.09    AIC: 113.5

The write-up:

Use the code in the previous section as well as your theoretical understanding of the variables to communicate how you chose and fit your starting model. Include:

  • what your error distribution assumption is and how you chose it

  • what your shape assumption is and how you chose it

  • what type of model you chose for your starting model (e.g. GLM)

For example,

I built a model to test my hypothesis that varibility in adult moose growth rate (delMass, kg per year) is explained by latitude (Lat, ˚N), minimum winter temperature (WTemp, ˚C) and sex (Sex, male or female). My error distribution assumption was normal (gaussian) as delMass is continuous and could be positive or negative. My shape assumption regarding the relationship between delMass and each covariate was linear. This was chosen after inspecting plots of delMass and both WTemp and Lat (Figures 1 & 2). Plots were made using the ggplot2 package (Wickham 2016)32 In addition, a linear shape assumption was chosen to describe the relationship between Sex and delMass as Sex is a categorical variable.

Based on these assumption I tested my hypothesis by fitting a generalized linear model with normal error distribution assumption to my data. All model fitting and analysis was done in R (R Core Team 2022).33

6.6 Assessing fit

The task:

Here, you will describe how you assessed your model to ensure that it can be used to test your hypothesis. Remember that the assumptions you made above (error distribution & shape assumptions) were a “best guess” but it is only after the model is fit that we can confirm that it is useful. In this section you will:

  • Consider covariate collinearity

  • Consider observation independence

  • Consider your error distribution assumption

  • Consider your shape assumption

and report your results.

The code:

Code supporting your tasks in this section will help you assess your model for violations that may make it invalid for testing your hypothesis, e.g.:

# Consider covariate collinearity 
## Fit a model without interactions
startMod.noInt <- glm(formula = delMass ~ Lat + WTemp + Sex, # hypothesis
              data = myDat, # data
              family = gaussian(link = "identity")) # error distribution assumption

startMod.noInt

Call:  glm(formula = delMass ~ Lat + WTemp + Sex, family = gaussian(link = "identity"), 
    data = myDat)

Coefficients:
(Intercept)          Lat        WTemp     SexMales  
    -13.013        0.879       -0.105       12.000  

Degrees of Freedom: 23 Total (i.e. Null);  20 Residual
Null Deviance:      1200 
Residual Deviance: 99.43    AIC: 112.2
## Calculate Variance Inflation Factors (VIFs)
#library(car) # load the car package
#vif(startMod.noInt) # estimate VIFs

## there are VIF values higher than the threshold value (in this case VIFs > 3).  Let's remove WTemp and recalculate the VIFs:

startMod.noInt.noWTemp <- glm(formula = delMass ~ Lat + Sex, # hypothesis
              data = myDat, # data
              family = gaussian(link = "identity")) # error distribution assumption

startMod.noInt.noWTemp

Call:  glm(formula = delMass ~ Lat + Sex, family = gaussian(link = "identity"), 
    data = myDat)

Coefficients:
(Intercept)          Lat     SexMales  
   -19.5975       0.9963      12.0000  

Degrees of Freedom: 23 Total (i.e. Null);  21 Residual
Null Deviance:      1200 
Residual Deviance: 99.58    AIC: 110.3
## Recalculating the VIFs

#vif(startMod.noInt.noWTemp) # estimate VIFs

## there are no longer any problems with covariate collinearity, as all VIFs < 3

## Refitting your starting model with the interactions back in:
startMod <- glm(formula = delMass ~ Lat + Sex + Lat:Sex, # hypothesis
              data = myDat, # data
              family = gaussian(link = "identity")) # error distribution assumption



# Consider observation independence
## You wonder if Source could be violating your observation independence assumption.  Check to see with a plot of your residuals vs. Source

library(DHARMa) # load package

simulationOutput <- simulateResiduals(fittedModel = startMod, n = 250) # simulate data from our model n times and calculate residuals

myDat$Resid <- simulationOutput$scaledResiduals # add residuals to data frame

ggplot()+ # start ggplot
  geom_violin(data = myDat,
             mapping = aes(x = Source, y = Resid))+ # add observations as a violin
  geom_point(data = myDat,
             mapping = aes(x = Source, y = Resid))+ # add observations as points
  xlab("Source university")+ # y-axis label
  ylab("Scaled residual")+ # x-axis label
  labs(caption = "Figure 3: A comparison of model residuals vs. university from which the data were sourced")+ # figure caption
  theme_bw()+ # change theme of plot
  theme(plot.caption = element_text(hjust=0)) # move figure legend (caption) to left alignment. Use hjust = 0.5 to align in the center.

  ## The pattern is uniform.  You assume no problem with observation dependence.

# Consider your error distribution assumption by inspecting residuals plotted vs. the fitted values
plot(simulationOutput, asFactor=FALSE) # compare simulated data to our observations

# Consider your shape assumption by inspecting residuals plotted vs. each covariate
plot(simulationOutput, # compare simulated data to 
     form=myDat$Sex, # our observations
     asFactor=TRUE) # whether the variable plotted is a factor

plot(simulationOutput, # compare simulated data to 
     form=myDat$Lat, # our observations
     asFactor=FALSE) # whether the variable plotted is a factor

The write-up:

Use the code in the previous section to comment on your model’s validity in testing your research hypothesis. Include:

  • how you determined if there were problems with covariate collinearity and any actions you took if you detected a problem

  • how you determined if there were problems with observation dependence and any actions you took if you detected a problem

  • how you determined if your error distribution assumption was valid and any actions you took to address problems

  • how you determined if your shape assumption was valid and any actions you took to address problems

For example,

I tested if collinearity among my covariates was making my model coefficients uncertain by estimating variance inflation factors (VIFs) with the car package (Fox & Weisberg 2019)34. Initial VIFs for Lat and WTemp were > 23 indicating a high level of covariate collinearity. WTemp was removed from the model and the VIFs were re-estimated. VIFs in the new model were both 1, and it was concluded that there was no further issue with covariate collinearity. The new starting model fits the hypothesis:

delMass ~ Latitude + Sex + Latitude:Sex

I tested my assumption of observation independence by determining if the observations were grouped by Source (the source university for the data35). I estimated my model residuals using the DHARMA package (Hartig 2022)36 and plotted my residuals vs. Source. I concluded my data were not dependent on one another (grouped by) Source as the residuals were evenly distributed across the three source universities (see figure 3).

I assessed my error distribution assumption by inspected my residuals. Observed residuals were similar to those expected given my normal error distribution assumption. The Kolmogorov-Smirnov test comparing the observed to the expected distribution was not significant (P = 0.96). The dispersion and outlier tests were also not significant (P = 0.74 and P = 0.99 respectively). A plot of the residuals vs. fitted values showed a uniform cloud. From these results, I concluded that my error distribution assumption was appropriate.

I assessed my shape assumption by inspected my residuals vs. each covariate. My residuals were evenly distributed with Latitude, indicating a linear shape assumption for Latitude was appropriate. The linear shape assumption for Sex was necessary as Sex is a categorical variable. A plot of the residuals vs. Sex showed residuals were evenly distributed across the two sexes. From these results, I concluded that my linear shape assumptions were appropriate.

Given these results, I determined that the new starting model (delMass ~ Latitude + Sex + Latitude:Sex) can be used to test my hypothesis.

6.7 Hypothesis test

The task:

Here, you will use your model to test your hypothesis. For the purposes of this course, you will do this with the model selection method we practiced in class37. In this section you will:

  • fit and compare models representing all possible combinations of the covariates in your starting model

  • use the results to choose best-specified model(s)

  • report your results.

The code:

Code supporting your tasks in this section will help you test and rank your models, e.g.:

library(MuMIn) # load package

options(na.action = "na.fail") # needed for dredge() function to prevent illegal model comparisons

(dredgeOut<-dredge(startMod, extra = "R^2")) # fit and compare a model set representing all possible predictor combinations
Global model call: glm(formula = delMass ~ Lat + Sex + Lat:Sex, family = gaussian(link = "identity"), 
    data = myDat)
---
Model selection table 
    (Int)    Lat Sex Lat:Sex    R^2 df  logLik  AICc delta weight
8  -1.596 0.7041   +       + 0.9340  5 -48.387 110.1  0.00  0.756
4 -19.600 0.9963   +         0.9170  4 -51.130 112.4  2.26  0.244
3  41.780          +         0.7200  3 -65.726 138.7 28.54  0.000
2 -13.600 0.9963             0.1971  3 -78.366 163.9 53.82  0.000
1  47.780                    0.0000  2 -80.999 166.6 56.46  0.000
Models ranked by AICc(x) 
## Our best-specified model is delMass ~ Lat + Sex + Lat:Sex (model #8 in the dredge() output).

bestMod<-(eval(attr(dredgeOut, "model.calls")$`8`)) # extract model #8 from dredge table

# You can make a pretty table to use in your write-up
library(gt) # load gt package

myTable <- gt(dredgeOut) # make a pretty table
myTable <- fmt_number(myTable, # to format the numbers in my table
                    columns = everything(), # which columns to format
                    decimals = 2) # round to 2 decimal places
myTable <- tab_caption(myTable, 
                       caption = "Table 1: Model selection table for hypothesis testing. Each row is a model fit to my data. Covariates included in each model are indicated by a number (for continuous covariates and intercept) or '+' (for categorical covariates). R^2 is the likelihood ratio R-squared, df indicates number of model parameters, logLik is the model Log-likelihood, AICc is the corrected Akaike Information Criteria, delta is the change in AICc between the model and that of the lowest AICc, and weight is the Akaike weights.")

myTable
Table 1: Model selection table for hypothesis testing. Each row is a model fit to my data. Covariates included in each model are indicated by a number (for continuous covariates and intercept) or '+' (for categorical covariates). R^2 is the likelihood ratio R-squared, df indicates number of model parameters, logLik is the model Log-likelihood, AICc is the corrected Akaike Information Criteria, delta is the change in AICc between the model and that of the lowest AICc, and weight is the Akaike weights.
(Intercept) Lat Sex Lat:Sex R^2 df logLik AICc delta weight
−1.60 0.70 + + 0.93 5.00 −48.39 110.11 0.00 0.76
−19.60 1.00 + NA 0.92 4.00 −51.13 112.36 2.26 0.24
41.78 NA + NA 0.72 3.00 −65.73 138.65 28.54 0.00
−13.60 1.00 NA NA 0.20 3.00 −78.37 163.93 53.82 0.00
47.78 NA NA NA 0.00 2.00 −81.00 166.57 56.46 0.00
## Could also be written as:
 # dredgeOut %>%
 #   gt() %>% # make a pretty table
 #   fmt_number(columns = everything(), # which columns to format
 #                     decimals = 2) %>% # round to 2 decimal places
 #   tab_caption(caption = "Table 1: Model selection table for hypothesis testing. Each row is a model fit to my data. Covariates included in each model are indicated by a number (for continuous covariates and intercept) or '+' (for categorical covariates). R^2 is the likelihood ratio R-squared, df indicates number of model parameters, logLik is the model Log-likelihood, AICc is the corrected Akaike Information Criteria, delta is the change in AICc between the model and that of the lowest AICc, and weight is the Akaike weights.")

The write-up:

Use the code in the previous section to explain how you chose your best-specified model. Include:

  • the method you are using to test your hypothesis

  • how you fit and rank your candidate model set

  • how you made your decision regarding your best-specified model

For example,

I used model selection to test my hypothesis that varibility in adult moose growth rate is explained by latitude, sex and the interaction between the two. I used the dredge() function from the MuMIn package (Bartón 2022) to fit and rank models representing all possible covariate combinations. Models were ranked by corrected Akaike Information Criteria38 which [add a brief description of what AICc is].39

The model with the lowest AICc was chosen as the best-specified model, with models within ∆AICc = 2 of the lowest AICc model being considered equally best-specified.

The best-specified model was the full model:

delMass ~ Lat + Sex + Lat:Sex

with an Akaike weight (normalized relative likelihood) of 0.756. The next best model had an AIC of 2.26 more than the top model and an Akaike weight of 0.244.

6.8 Reporting

The task:

In this section, you will report your hypothesis testing results. Specifically,

  • Report your best-specified model(s)

  • Report your modelled effects: This includes both reporting your coefficients and plotting the effects. I’ve moved visualization here to be more consistent with how we report biological results. A few things to note here:

    • reporting the coefficients includes:

      • reporting your coefficient estimates with uncertainty: make sure you report appropriate significant digits and units.

      • reporting evidence that these coefficient estimates are different than zero: the fact that your covariate is in your best-specified model means that there is evidence of an effect of the covariate on the response (i.e. that the coefficient associated with the covariate is different than zero), BUT if you have a categorical variable with multiple levels, you need to test to see which coefficients for each level are different than zero.

      • reporting evidence that these coefficient estimates are different than each other across the levels of your categorical covariate (if you have one): this is the multiple comparison testing we practiced in class. Remember that if you have both a continuous covariate and a categorical covariate, you need to test the coefficients (slopes) for the continuous covariate first and only if the slopes are similar across levels of the categorical covariate can you test if the intercepts are significantly different. This is due to the fact that we can get different intercept estimates just due to slopes not being equal.

      • In your visualization, remember to include your modelled fit, uncertainty and your observations in your plot)

  • Report how well your model explains your response (explained deviance)

  • Report the relative importance of each term in your model

The code:

Code supporting your tasks in this section will help you communicate what your results tell you about your hypothesis, e.g.:

# Reporting the modelled effects of your covariates on your response: 

## ** your coefficient estimates with uncertainty
### Find the coefficients for Latitude.  As Sex is categorical, we do this for emmeans() and the reported coefficients are predicted mean response at each level of Sex. Note that emmeans() can report coefficients on the RESPONSE scale.
library(emmeans) # load emmeans package
emmeans(object = bestMod, # your model
        specs = ~ Lat + Sex + Lat:Sex, # your covariates
       # at = list(Lat = 0), # to get coefficient predictions for Sex when Lat = 0
        type = "response") # report coefficients on the response scale
  Lat Sex     emmean    SE df lower.CL upper.CL
 61.6 Females   41.8 0.575 20     40.6       43
 61.6 Males     53.8 0.575 20     52.6       55

Confidence level used: 0.95 
### Find the coefficients for Latitude.  As Latitude is continuous, we do this with emtrends() and the reported coefficients are the slope describing the change in predicted mean response with a unit change in Latitude.  Note that emtrends() reports coefficients on the LINK scale.  You need to convert this to the response scale.  In our case, the error distribution assumption is normal so nothing needs to be done to convert the coefficients.
emtrends(object = bestMod, # your model
        specs = ~ Lat + Sex + Lat:Sex, # your covariates
        var = "Lat") # name of your continuous cvoariate
  Lat Sex     Lat.trend    SE df lower.CL upper.CL
 61.6 Females     0.704 0.182 20    0.324     1.08
 61.6 Males       1.289 0.182 20    0.908     1.67

Confidence level used: 0.95 
## ** evidence that these coefficient estimates are different than zero & ** evidence that these coefficient estimates are different than each other across the levels of your categorical covariate (if you have one)

### Are the coefficients (slopes) of the latitude effect on growth different?
forLatCompare <- emtrends(object = bestMod, # your model
                  specs = pairwise ~ Lat + Sex + Lat:Sex, # your covariates
                  var = "Lat") # name of your continuous cvoariate

test(forLatCompare) # to get results of slope tests and comparisons
$emtrends
  Lat Sex     Lat.trend    SE df t.ratio p.value
 61.6 Females     0.704 0.182 20   3.861  0.0010
 61.6 Males       1.289 0.182 20   7.065  <.0001


$contrasts
 contrast                        estimate    SE df t.ratio p.value
 Lat61.6 Females - Lat61.6 Males   -0.584 0.258 20  -2.266  0.0347
plot(forLatCompare, # plot coefficients with error
     comparison = TRUE, # include arrows indicating thresholds for differences
     type = "response") # plot on the response scale

### Are the coefficients (intercepts) of Sex factor levels different than zero?
forSexCompare <- emmeans(object = bestMod, # your model
                         specs = pairwise ~ Lat + Sex + Lat:Sex, # your covariates
                         type = "response") # name of your continuous cvoariate


test(forSexCompare)
$emmeans
  Lat Sex     emmean    SE df t.ratio p.value
 61.6 Females   41.8 0.575 20  72.704  <.0001
 61.6 Males     53.8 0.575 20  93.588  <.0001


$contrasts
 contrast                        estimate    SE df t.ratio p.value
 Lat61.6 Females - Lat61.6 Males      -12 0.813 20 -14.768  <.0001
## Plot your modelled effects:
# Set up your covariates for the visualized fit
forLat<-seq(from = min(myDat$Lat), 
             to = max(myDat$Lat), 
             by = 1) # get a range of latitudes for making predictions
forSex<-unique(myDat$Sex) # get every level of my Sex covariate
forVis<-expand.grid(Lat=forLat, Sex=forSex) # create a data frame with all combinations of covariates

# Get your model fit estimates at each value of your covariates
modFit<-predict(bestMod, # the model
                newdata = forVis, # the covariate values
                type = "response", # make the predictions on the response scale
                se.fit = TRUE) # include uncertainty estimate

forVis$Fit<-modFit$fit # add your fit to the data frame
forVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frame
forVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame

ggplot()+
  geom_point(data = myDat, # data
             mapping = aes(x = Lat, y = delMass, col = Sex))+ # add data to your plot
  geom_ribbon(data = forVis, 
              mapping = aes(x = Lat, ymin = Lower, ymax = Upper, fill = Sex), alpha = 0.5)+ # add the uncertainty to your plot
  geom_line(data = forVis, 
            mapping = aes(x = Lat, y = Fit, col = Sex))+ # add the model fit to your plot
  ylab("Change in body mass (kg per year)")+ # change y-axis label
  xlab("Latitude (˚N)")+ # change x-axis label
  labs(caption = "Figure 4: Modelled effects (lines, mean +/- standard error) of latitude and sex (colours) on growth rate \nalong with observations (points)")+
  theme_bw()+ # change ggplot theme
  theme(plot.caption = element_text(hjust=0)) # move figure legend (caption) to left alignment. Use hjust = 0.5 to align in the center.

# * Report how well your model explains your response (explained deviance)

r.squaredLR(bestMod) # estimates the likelihood ratio R^2.  Could also estimate a traditional R^2 with 1-summary(bestMod)$deviance/summary(bestMod)$null.deviance here as the error distribution assumption is normal and the shape assumption is linear, but the likelihood ratio R^2 function is generally applicable to many error distribution assumptions and equivalent to the traditional R^2 when the error distribution assumption is normal and the shape assumption is linear.
[1] 0.9339728
attr(,"adj.r.squared")
[1] 0.9350677
# * report how important each model term (covariate fixed effects and any interactions) is in explaining the deviation in your response. 
sw(dredgeOut)
                     Sex  Lat  Lat:Sex
Sum of weights:      1.00 1.00 0.76   
N containing models:    3    3    1   

The write-up:

Use the results of the code above to report your hypothesis testing results. Include:

  • a description of the terms (covariate fixed effects and any interactions) that are in your best-specified model

  • compare these terms to your starting model - are all terms in your starting model in you best-specified model?

  • report your modelled effects by reporting your coefficients. As we will discuss this week in class, report:

    • your coefficient estimates with uncertainty

    • evidence that these coefficient estimates are different than zero

    • evidence that these coefficient estimates are different than each other across the levels of your categorical covariate (if you have one)

  • report your modelled effects by plotting. Remember to include your model predictions, uncertainty as well as your observations. Include units on your axes and a figure number and legend.

  • report your estimate of the % of deviance in your response that your model explains

  • consider what might be explaining the remaining (unexplained) deviation in your response

  • report how important each model term (covariate fixed effects and any interactions) is in explaining the deviation in your response.

For example,

My best-specified model tells me that the growth varies with latitude and sex, and that the effect of latitude on growth varies with sex. The terms in my best-specified model are the same as those in my starting model indicating that there is evidence that the main effects of latitude and sex as well as the interaction are explaining variability in moose growth.

Growth rate is higher for males than females: the predicted growth when latitude is 61.6˚N (the mean of the latitudinal range) is 41.8 ± 0.6 kg year-1 for females and 53.8 ± 0.6 kg year-1 for males.40 Both these coefficients are significantly different than zero (t-test, P < 0.0001 for males and females)^[note that we don’t need to add an adjustment for multiple comparisons here as there are only two levels in our categorical covariate).

The coefficient (slope) for the effect of a change of latitude on the growth rate is 0.70 ± 0.18 kg year-1 ˚N-1 for females and 1.29 ± 0.18 kg year-1 ˚N-1 for males. These slopes are significantly different from zero (t-test, P = 0.001 for females and P < 0.0001 for males) and significantly different from one another (t-test, P = 0.035). These coefficients show that growth rates increase with latitude for both sexes, but that the effect of latitude on growth is higher for male vs. female moose.

Figure 4 shows the modelled effects of latitude and sex on growth rate along with my observations.

Together the effects of latitude and sex explain 93% of the deviance in growth rate (Likelihood ratio R2).41

The remaining unexplained deviance may be due to other factors affecting growth rate such as winter harshness, food availability, etc. Note that I was not able to include minimum winter temperature in my hypothesis test due to high collinearity between latitude and minimum winter temperature. Therefore, I do not know if the measured latitudinal effect on growth might mechanistically be due to winter harshness (here estimated as minimum winter temperature). This could be tested if I was able to expand my data set to include sites where latitude and minimum winter temperature were less correlated.

Based on the sum of Akaike weights, latitude and sex are equally important in explaining the deviance in growth rate as they appear in all models with Akaike weights > 0. The interaction between latitude and sex is slightly less important appearing only in the best-specified model with an Akaike weight of 0.76.

6.9 Predicting

The task:

Here you may use your model to make predictions (remembering prediction limits). For the assignments and exam, I’ll explicitly ask you to make a prediction if we want to see one.

If you want to make a prediction (e.g. based on your best-specified model, what is your predicted mean annual growth of a female moose living in Sicily, Italy?), in this section you want to report your prediction results and give any limitations to the prediction.

The code:

Code supporting your tasks in this section will help you communicate what prediction you made and the results, e.g.:

# Based on your best-specified model, what is your predicted mean annual growth of an adult female moose living in Sicily, Italy?

## Looking on the internet, we see that Sicily, IT is at 37.6˚N.

predict(bestMod, # our model
        newdata = data.frame(Sex = "Females", Lat = 37.6), # values of the covariates at which to make the prediction
        se.fit=TRUE, # include an estimate of error around the prediction
        type="response") # make sure the prediction is on the response scale
$fit
       1 
24.87708 

$se.fit
[1] 4.414466

$residual.scale
[1] 1.990439

The write-up:

Use the code in the previous section to present the results of your prediction. Include

  • the prediction and an estimate of uncertainty

  • any perceived prediction limits.

For example,

I used my best-specified model to estimate the mean annual growth rate of adult female moose living in Sicily, Italy (~ 37.6˚N) as 24.9 +/- 4.4 kg per year (mean +/- standard error). While this is an estimate consistent with my model, it likely is unrealistic as the distribution of moose does not currently extend as far south as Sicily. It is likely there are limitations to their dispersal to or survival in this area and my model may not be valid for latitudes so far south.

- Reporting your best-specified model   - Reporting your modelled patterns      - Plotting model fit, uncertainty and observations [[Visualizing]]      - Reporting effect sizes        - Report the coefficients - [[what are model coefficients]]         - What are the coefficient estimates (with uncertainty)?            - Is each coefficient significantly different from zero?            - Are the coefficients significantly different than one another?    - Reporting how well your model explains your response (explained variation)    - Report the relative importance of each term in your model

cross validation

visreg vs. ggeffects

three ways to report what your model is saying about your hypothesis
- visualize
- report coefficients
- give an example (response predicted under conditions)


- model coefficients are model effects
    - the effect on Resp with a change in Cov
- Continuous
    - coefficient is the change in Resp with a unit change in Cov
    - what does it mean if coef is positive
    - what does it mean if coef is negative
    - what is error on the coef?
    - what does it mean when a coef is no different than 0? - why that is the science
    - how to test if two coefficients are similar
- Categorical
    - coefficient is the change in Resp with a change from one level of Cat to another
    - what does it mean if coef is positive
    - what does it mean if coef is negative
    - what is error on the coef?
    - what does it mean when a coef is no different than 0? - why that is the science
    - how to test if two coefficients are similar
- Show linear model equation and how it looks for Cat, Cont with coefficients
- Describe generically and then define for each link function

report effect size and precision, not p-value 


Emphasise effect sizes to replace statistical significance with ecological relevance
research question should follow the FINER criteria
refinements of your hypothesis after you collect the data will result in a hypothesis that reflects the sample, not the population
exploratory or hypothesis generating analysis
use representative (unbiased) sampling
to report results of hypothesis testing, report effect sizes, confidence intervals, interpretation of effect and confidence intervals, and can include p-values
p-values are the probability that we would obtain a result equal to or more extreme than that observed given a null hypothesis that is true
p-values tell us more about how compatible the data are with the model vs. evidence of underlying effects

7 A statistical modelling example

It will help our discussions to have an example of a best-specified model. For our example, let’s assume that your hypothesis is that your response variable Resp is explained by a categorical covariate (Cat), a continuous covariate (Cont) as well as the interaction between the two. You can communicate this as:

Resp ~ Cat + Cont + Cat:Cont

Let’s assume you tested this hypothesis by fitting a model with a Gamma error distribution assumption (to reflect the nature of your response variable) and linear shape assumption (to reflect the relationship among each covariate and your response) using a GLM:

str(myDat) # a look at our data
'data.frame':   50 obs. of  3 variables:
 $ Cont: int  50 64 65 54 52 39 53 41 44 34 ...
 $ Cat : Factor w/ 2 levels "Control","Treatment": 2 1 1 2 1 2 1 2 2 2 ...
 $ Resp: num  27.7 17.2 20.1 39.4 15.4 ...
startMod<-glm(formula = Resp ~ Cont + Cat + Cont:Cat, # hypothesis
               data = myDat, # data
               family = Gamma(link="inverse")) # error distribution assumption

You considered possible covariate collinearity42 and observation independence43 and how it might be affecting your starting model. You then assessed your starting model to make sure the error distribution and shape assumptions were appropriate using the DHARMa method:

library(DHARMa)
simulationOutput <- simulateResiduals(fittedModel = startMod, n = 250) # simulate data from our model n times
 
plot(simulationOutput, asFactor=FALSE) # compare simulated data to our observations

plot(simulationOutput, # compare simulated data to 
     form=myDat$Cont, # our observations
     asFactor=FALSE) # whether the variable plotted is a factor

plot(simulationOutput, # compare simulated data to 
     form=myDat$Cat, # our observations
     asFactor=TRUE) # whether the variable plotted is a factor

All looks fine - the distribution of the expected and observed residuals are similar, and there is no pattern in the residuals when plotted with fitted values or any of the covariates.

You then used your starting model to test your hypothesis by model selection using the dredge() function in the MuMIn package:

library(MuMIn)

options(na.action = "na.fail") # needed for dredge() function to prevent illegal model comparisons

(dredgeOut<-dredge(startMod, extra = "R^2")) # fit and compare a model set representing all possible predictor combinations
Global model call: glm(formula = Resp ~ Cont + Cat + Cont:Cat, family = Gamma(link = "inverse"), 
    data = myDat)
---
Model selection table 
    (Int) Cat        Cnt Cat:Cnt    R^2 df   logLik  AICc delta weight
8 0.07407   + -0.0003389       + 0.8703  5 -123.363 258.1  0.00  0.815
4 0.08742   + -0.0006263         0.8554  4 -126.083 261.1  2.97  0.185
2 0.05874   +                    0.6822  3 -145.763 298.0 39.96  0.000
3 0.08423     -0.0008670         0.3726  3 -162.769 332.1 73.97  0.000
1 0.04133                        0.0000  2 -174.425 353.1 95.02  0.000
Models ranked by AICc(x) 

You find two models have something to say about your hypothesis but the top model (model #8) is greater than 2 AICc lower than the second model (model #4). You choose your best-specified model as:

bestMod<-(eval(attr(dredgeOut, "model.calls")$`8`)) # extract model #8 from dredge table

bestMod

Call:  glm(formula = Resp ~ Cat + Cont + Cat:Cont + 1, family = Gamma(link = "inverse"), 
    data = myDat)

Coefficients:
      (Intercept)       CatTreatment               Cont  CatTreatment:Cont  
        0.0740663         -0.0025851         -0.0003389         -0.0004083  

Degrees of Freedom: 49 Total (i.e. Null);  46 Residual
Null Deviance:      5.916 
Residual Deviance: 0.7805   AIC: 256.7

Here we’ll use this example best-specified model to visualize and report what it is saying about your hypothesis. We’ll cover other examples in class as well.

8 Visualizing

Visualizing your model means visualizing what your model is saying about your research hypothesis. We have practiced visualizing your fitted statistical model but now let’s focus on visualizing to communicate the effects described in your model.

Here, you want to visualize how each covariate is affecting the response by drawing how fitted values (the estimates of your response made by the model) change with your covariates, along with a measure of uncertainty around those fitted values. You also want to include your data (actual observations of your response and covariates) so you can start to communicate how well your model captures your hypothesis, and what amount of deviance in your response is explained by your model.

There are a lot of different R packages that make it easy to quickly visualize your model. We’ll focus on two methods that will allow you to make quick visualizations of your model and/or customize the figures as you would like.

8.1 Visualizing with the visreg package

The first method uses the visreg package and the visreg() function which we’ve practiced before. Using your example bestMod, here are a few visualizations using the visreg package:

Visualizing with a continuous covariate on the x-axis: Notice how you can include the gg = TRUE argument to plot this as a ggplot type figure. This allows you to add your data onto the visualization of your model.

library(visreg) # load visreg package
library(ggplot2) # load ggplot2

visreg(bestMod, # model to visualize
       scale = "response", # plot on the scale of the response
       xvar = "Cont", # covariate on x-axis
       by = "Cat", # covariate plotted as colour
       #breaks = 3, # if you want to control how the colour covariate is plotted
       #cond = , # if you want to include a 4th covariate
       overlay = TRUE, # to plot as overlay or panels 
       rug = FALSE, # to include a rug
       gg = TRUE)+ # to plot as a ggplot
  geom_point(data = myDat, # data
             mapping = aes(x = Cont, y = Resp, col = Cat))+ # add data to your plot
  ylim(0, 60)+ # adjust the y-axis units
  ylab("Response, (units)")+ # change y-axis label
  xlab("Cont, (units)")+ # change x-axis label
  theme_bw() # change ggplot theme

Visualizing with a categorical covariate on the x-axis: You can also plot the same model with the categorical covariate on the x-axis. Notice how the continuous covariate is represented by different levels in the colours on the plot. Here you’ve asked for three levels with breaks = 3.

visreg(bestMod, # model to visualize
       scale = "response", # plot on the scale of the response
       xvar = "Cat", # covariate on x-axis
       by = "Cont", # covariate plotted as colour
       breaks = 3, # if you want to control how the colour covariate is plotted
       #cond = , # if you want to include a 4th covariate
       overlay = TRUE, # to plot as overlay or panels 
       rug = FALSE, # to include a rug
       gg = TRUE)+ # to plot as a ggplot
  ylab("Response, (units)")+ # change y-axis label
  xlab("Cat, (units)")+ # change x-axis label
  theme_bw() # change ggplot theme

You can also specify at which levels the breaks should occur with the breaks = ... argument. For example, you can ask visreg to plot the modelled effects when Cont = 40 and Cont = 50:

visreg(bestMod, # model to visualize
       scale = "response", # plot on the scale of the response
       xvar = "Cat", # covariate on x-axis
       by = "Cont", # covariate plotted as colour
       breaks = c(40,50), # if you want to control how the colour covariate is plotted
       #cond = , # if you want to include a 4th covariate
       overlay = TRUE, # to plot as overlay or panels 
       rug = FALSE, # to include a rug
       gg = TRUE)+ # to plot as a ggplot
  ylab("Response, (units)")+ # change y-axis label
  xlab("Cat, (units)")+ # change x-axis label
  theme_bw() # change ggplot theme

Unfortunately due to limitations of the visreg package, you can not easily add our data onto plots where the x-axis is a categorical covariate. But that’s ok, because there are other options…

8.2 Visualizing “by hand”:

We will also practice how to visualize your model “by hand”. Here, “by hand” is a bit of a silly description as R will be doing the work for you. What I mean by “by hand” is that you will be building the plot yourselves by querying your model object.

Visualizing with a continuous covariate on the x-axis:. To plot by hand, you will first make a data frame containing the value of your covariates at which you want to plot effects on the response:

# Set up your covariates for the visualized fit
forCont<-seq(from = min(myDat$Cont), 
             to = max(myDat$Cont), 
             by = 0.1) # e.g. a sequence of Cont values
forCat<-unique(myDat$Cat) # every value of your categorical covariate

# create a data frame with your covariates
forVis<-expand.grid(Cont=forCont, Cat=forCat) # expand.grid() function makes sure you have all combinations of covariates

Next, you will use the predict() function44 to get the model estimates of your response variable at those values of your covariates:

# Get your model fit estimates at each value of your covariates
modFit<-predict(bestMod, # the model
                newdata = forVis, # the covariate values
                type = "response", # make the predictions on the response scale
                se.fit = TRUE) # include uncertainty estimate

forVis$Fit<-modFit$fit # add your fit to the data frame
forVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frame
forVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame

Finally, you will use these model estimates to make your plot:

ggplot()+
  geom_point(data = myDat, # data
             mapping = aes(x = Cont, y = Resp, col = Cat))+ # add data to your plot
  geom_ribbon(data = forVis, 
              mapping = aes(x = Cont, ymin = Lower, ymax = Upper, fill = Cat), alpha = 0.5)+ # add the uncertainty to your plot
  geom_line(data = forVis, 
              mapping = aes(x = Cont, y = Fit, col = Cat))+ # add the model fit to your plot
  ylim(0, 60)+ # adjust the y-axis units
  ylab("Response, (units)")+ # change y-axis label
  xlab("Cont, (units)")+ # change x-axis label
  theme_bw() # change ggplot theme

Visualizing with a categorical covariate on the x-axis:. The procedure for visualizing when you have a categorical covariate on the x-axis is similar: 1) choose the values of your covariates at which to make a prediction, 2) use predict() to use your model to estimate your response variable at those values of your covariate, and 3) use the model estimates to plot your model fit:

# Set up your covariates for the visualized fit
forCont<-c(40, 50 ,60) # e.g. Cont values
forCat<-unique(myDat$Cat) # every value of your categorical covariate

# create a data frame with your covariates
forVis<-expand.grid(Cont=forCont, Cat=forCat) # expand.grid() function makes sure you have all combinations of covariates

# Get your model fit estimates at each value of your covariates
modFit<-predict(bestMod, # the model
                newdata = forVis, # the covariate values
                type = "response", # make the predictions on the response scale
                se.fit = TRUE) # include uncertainty estimate

forVis$Fit<-modFit$fit # add your fit to the data frame
forVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frame
forVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame


ggplot()+
  geom_jitter(data = myDat, # data
             mapping = aes(x = Cat, y = Resp))+ # add data to your plot
  geom_errorbar(data = forVis, 
              mapping = aes(x = Cat, y = Fit, col = Cont, ymin = Lower, ymax = Upper))+ # add the uncertainty to your plot
  ylim(0, 60)+ # adjust the y-axis units
  ylab("Response, (units)")+ # change y-axis label
  xlab("Cat, (units)")+ # change x-axis label
  theme_bw() # change ggplot theme

8.3 Plotting in 3D

You might notice the plots above are communicating 3-dimensions (one response + two covariates) in a 2-dimensional plot. There are other ways of making 3-dimensional plots in R, e.g. with the visreg package using the visreg2d() function in the visreg package:

visreg2d(fit = bestMod, # your model
         xvar = "Cont", # what to plot on the x-axis
         yvar = "Cat", # what to plot on the y-axis
         scale = "response") # make sure fitted values (colours) are on the scale of the response

or “by hand” using the geom_raster() function in the ggplot2 package:

# Set up your covariates for the visualized fit
forCont<-seq(from = min(myDat$Cont), 
             to = max(myDat$Cont), 
             by = 0.1) # e.g. a sequence of Cont values
forCat<-unique(myDat$Cat) # every value of your categorical covariate

# create a data frame with your covariates
forVis<-expand.grid(Cont=forCont, Cat=forCat) # expand.grid() function makes sure you have all combinations of covariates


# Get your model fit estimates at each value of your covariates
modFit<-predict(bestMod, # the model
                newdata = forVis, # the covariate values
                type = "response", # make the predictions on the response scale
                se.fit = TRUE) # include uncertainty estimate

forVis$Fit<-modFit$fit # add your fit to the data frame
forVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frame
forVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame

# create your plot
ggplot() + # start your ggplot
  geom_raster(data = forVis, aes(x = Cont, y = Cat, fill = Fit))+ # add the 3 dimensions as a raster
  geom_point(data = myDat, aes(x = Cont, y = Cat, colour = Resp)) # add your data

EXTRA: As you can see, these plots represent the 3rd dimension by using colour. We can also make actual 3 dimensional plots in R with the plotly package45. These plots are interactive which makes them more useful than static 3d plots. Use your mouse to move the plot around!

library(plotly) # load the plotly package

plot_ly(data = forVis, # the data with your model predictions (made above)
        x = ~Cont, # what to plot on the x axis
        y = ~Cat, # what to plot on the y axis
        z = ~Fit, # what to plot on the z axis
        type="scatter3d", # type of plot
        mode="markers") %>% # type of plot
  add_markers(data = myDat, x = ~Cont, y = ~Cat, z = ~Resp) # add your data

9 Reporting

To report your model, you want to communicate a number of things:
* Reporting your best-specified model(s). * Reporting the modelled effects of your covariates on your response.
* Reporting how well your model explains your response. * Reporting the importance of each covariate in explaining your response.

9.1 Reporting your best-specified model(s)

Reporting your best-specified model means reporting the terms - the covariates and any interactions - that are in your model. It’s good practice to present the model along with the results from the model selection. In this way, you can include multiple best-specified models if there is evidence that more than one might be useful.

For example, with your results from model selection being:

Global model call: glm(formula = Resp ~ Cont + Cat + Cont:Cat, family = Gamma(link = "inverse"), 
    data = myDat)
---
Model selection table 
    (Int) Cat        Cnt Cat:Cnt    R^2 df   logLik  AICc delta weight
8 0.07407   + -0.0003389       + 0.8703  5 -123.363 258.1  0.00  0.815
4 0.08742   + -0.0006263         0.8554  4 -126.083 261.1  2.97  0.185
2 0.05874   +                    0.6822  3 -145.763 298.0 39.96  0.000
3 0.08423     -0.0008670         0.3726  3 -162.769 332.1 73.97  0.000
1 0.04133                        0.0000  2 -174.425 353.1 95.02  0.000
Models ranked by AICc(x) 

You can present the top two models as:

Depending on your hypothesis and results, you may want to present all models within ∆AICc < 2 of the best-specified model, or all models with any Akaike weight (as above), or simply all models. We’ll discuss ways of making these choices in class.

9.2 Reporting the modelled effects of your covariates on your response.

You want to report how the covariates are affecting your response. e.g. does a change in the covariate lead to an increase in the response? or a decrease? or is the effect non-linear (for continuous covariates)?

The easiest way to do this is by visualization, as you saw above, but you also want to report the coefficients with their uncertainty:

summary(bestMod)$coefficients
                       Estimate   Std. Error    t value     Pr(>|t|)
(Intercept)        0.0740662767 0.0070213172 10.5487724 7.241364e-14
CatTreatment      -0.0025850597 0.0087148893 -0.2966257 7.680875e-01
Cont              -0.0003388666 0.0001490856 -2.2729675 2.774499e-02
CatTreatment:Cont -0.0004082566 0.0001772830 -2.3028518 2.586062e-02

Note that for categorical covariates with more than two levels (e.g. an example with a covariate consisting of Species A, B and C), you’ll want to take an extra step to determine if the effect of each level of the covariate is similar (e.g. is the effect of Species A on the response the same as the effect on Species B?). We’ll cover how to do this in Week 12.

9.3 Reporting how well your model explains your response

If you will recall, your whole motivation for pursuing statistical modelling was to explain variation in your response. Thus, it is important that you quantify how much variation in your response you are able to explain by your model.

Note that we will discuss this as “explained deviance” rather than “explained variation”. This is because the term “variance” is associated with models where the error distribution assumption is normal, whereas deviance is a more universal term.

When you have a normal error distribution assumption and linear shape assumption, you can capture the amount of explained deviance simply by comparing the variance in your response (null variation) vs. the variance in your model residuals (residual variation) as the \(R^2\):

\(R^2 = 1 - \frac{residual variation}{null variation}\)

From this equation, you can see how, if your model is able to explain all the variation in the response, the residual variation will be zero, and \(R^2 = 1\). Alternatively, if the model explains no variation in the response the residual variation equals the null variation and \(R^2 = 0\).

For models with other error distribution and shape assumptions, you need another way of estimating the goodness of fit of your model. You can do this through a pseudo \(R^2\).

One useful pseudo \(R^2\) is called the Likelihood Ratio \(R^2\) or \(R^2_{LR}\). The \(R^2_{LR}\) compares the likelihood of your best-specified model to the likelihood of the null model:

\(R^2_{LR} = 1-exp(-\frac{2}{n}(log𝓛(model)-log𝓛(null)))\)

where \(n\) is the number of observations, \(log𝓛(model)\) is the log-likelihood of your model, and \(log𝓛(null)\) is the log-likelihood of the null model. The \(R^2_{LR}\) is the type of pseudo \(R^2\) that shows up in your dredge() output when you add extra = "R^2" to the dredge() call. You can calculate \(R^2_{LR}\) by hand, read it from our dredge() output, or estimate it using r.squaredLR() from the MuMIn package:

r.squaredLR(bestMod)
[1] 0.8702944
attr(,"adj.r.squared")
[1] 0.8711072

Note here that two values of \(R^2_{LR}\) are reported. The adjusted pseudo \(R^2\) given here under attr(,"adj.r.squared") has been scaled so that \(R^2_{LR}\) can reach a maximum of 1, to be equivalent to a regular \(R^2\).

One nice feature of the \(R^2_{LR}\) is that it is equivalent to the regular \(R^2\) when our model assumes a normal error distribution assumption and linear shape assumption, so you can use \(R^2_{LR}\) for any of the models we’re discussing in class.

  • check RR2 package

9.4 Reporting the importance of each covariate in explaining your response.

Finally, you want to report how relatively important each term is in your model in regards to explaining deviance in your response. This is also called “partitioning the explained deviance” to the covariates.

To get an estimate of how much deviance in your response one particular covariate explains, you may be tempted to compate the explained deviance (\(R^2\)) estimates of models fit to data with and without that particular covariate. Let’s try an example using your example statistical model:

dredgeOut
Global model call: glm(formula = Resp ~ Cont + Cat + Cont:Cat, family = Gamma(link = "inverse"), 
    data = myDat)
---
Model selection table 
    (Int) Cat        Cnt Cat:Cnt    R^2 df   logLik  AICc delta weight
8 0.07407   + -0.0003389       + 0.8703  5 -123.363 258.1  0.00  0.815
4 0.08742   + -0.0006263         0.8554  4 -126.083 261.1  2.97  0.185
2 0.05874   +                    0.6822  3 -145.763 298.0 39.96  0.000
3 0.08423     -0.0008670         0.3726  3 -162.769 332.1 73.97  0.000
1 0.04133                        0.0000  2 -174.425 353.1 95.02  0.000
Models ranked by AICc(x) 

If you want to get an estimate as to how much response deviance the Cont covariate explains, you might compare the pseudo \(R^2\) of a model with and without the Cont covariate, e.g. comparing model #4 (that includes Cont) and model #2 (that doesn’t include Cont):

# deviance explained by Cont:
R2.mod4 <- (dredgeOut$`R^2`[2]) # model #4 is the second row in dredgeOut
R2.mod2 <- (dredgeOut$`R^2`[3]) # model #2 is the third row in dredgeOut

R2.mod4 - R2.mod2 # find estimated contribution of Cont to explained deviance
[1] 0.1731384

But what if you chose to compare model #3 (that includes Cont) and model #1 (that doesn’t include Cont):

# deviance explained by Cont:
R2.mod3 <- (dredgeOut$`R^2`[4]) # model #3 is the fourth row in dredgeOut
R2.mod1 <- (dredgeOut$`R^2`[5]) # model #1 is the fifth row in dredgeOut

R2.mod3 - R2.mod1 # find estimated contribution of Cont to explained deviance
[1] 0.372649

Quite a different answer! Your estimates of the contribution of Cont to explaining the response deviation don’t agree because of collinearity among our covariates46. We’ll discuss this more in class.

One option for partitioning the explained deviance when you have collinearity among your covariates is hierarchical partitioning. Hierarchical partitioning estimates the average independent contribution of each covariate to the total explained variance by considering all models in the candidate model set47. This method is beyond the scope of the course but see the rdacca.hp package for an example of how to do this.

Another method for estimating the importance of each term (covariate or interaction) in your model is by again looking at the candidate model set ranking made by dredge(). Here you can measure the importance of a covariate by summing up the Akaike weights for any model that includes a particular covariate. The Akaike weight of a model compares the likelihood of the model scaled to the total likelihood of all the models in the candidate model set. The sum of Akaike weights for models including each covariate tells us how important the covariate is in explaining the deviance in your response. You can calculate the sum of Akaike weights with the sw() function in the MuMIn package:

sw(dredgeOut)
                     Cat  Cont Cat:Cont
Sum of weights:      1.00 1.00 0.82    
N containing models:    3    3    1    

Here we can see that Cat and Cont are equally important in explaining the deviance in Resp (they appear in all models that have any Akaike weight), while the interaction term between Cat and Cont is less important (only appearing in one model with Akaike weight, though this is the top model).

Check relaimpo

10 Predicting

Now that you have your best-specified model, you may want to use it to make a prediction, e.g. “what would I expect the value of my response to be at a different value of my covariate?” This may include estimating your response at a different place, time or environmental (covariate) condition.

10.1 First, consider the limits to making predictions

The first step when considering using your model to make predictions is to decide if this is appropriate. By using your model to make a prediction, you are assuming that the same processes that govern the relationship between your covariate(s) and response will hold in the new place, time or environmental condition. For example, you might have used a linear shape assumption that proved appropriate for modelling how growth varies due to temperature, but now you want to estimate growth at an even higher temperature than you tested. Is this appropriate? Would it be possible that growth starts to respond non-linearly with temperature?

10.2 How to use your model to make predictions

Once you are convinced that making a prediction with your model is useful, you can use the predict() function to predict the value of your response by handing it values of each covariate at which you would like a response prediction. For example:

predict(object = bestMod, # your model
        newdata = data.frame(Cont = 57, Cat = "Treatment"), # the values of the covariates at which to make the prediction
        type = "response", # to make the prediction on the response scale.
        se.fit = TRUE) # to include a measure of uncertainty around the prediction
$fit
       1 
34.60783 

$se.fit
       1 
1.060142 

$residual.scale
[1] 0.129131

So with Cont = 57 and Cat = Treatment, we would expect Resp to be 34.6 ± 1.1. This can be shown on our visualization here (prediction in black):

11 Reporting models with categorical covariates

In this week’s notes, we’ll work through three example best-specified models to report our model results, and in some cases make predictions.

Note that each example is testing a different hypothesis by fitting a different model to a different data-set. The three examples are not related to one another.

Note that I’m not going through all the steps that got us to these example best-specified models, but assume we went through the previous steps (Response, Covariates, Hypothesis, Starting model, Assessing Fit and Hypothesis Test) as we practiced in class to choose each best-specified model.

11.1 Example 1: Resp1 ~ Cat1

For example #1, assume your best-specified model shows that your response variable (Resp1) is explained by a categorical covariate (Cat1):

Resp1 ~ Cat1

Your best-specified model for example #1 is in an object called bestMod1:

bestMod1

Call:  glm(formula = Resp1 ~ Cat1 + 1, family = gaussian(link = "identity"), 
    data = myDat1)

Coefficients:
(Intercept)      Cat1Sp2      Cat1Sp3  
     22.500        4.194       -5.889  

Degrees of Freedom: 99 Total (i.e. Null);  97 Residual
Null Deviance:      8882 
Residual Deviance: 7037     AIC: 717.2

that was fit to data in myDat1:

summary(myDat1)
  Cat1        Resp1      
 Sp1:28   Min.   : 1.00  
 Sp2:36   1st Qu.:15.00  
 Sp3:36   Median :22.00  
          Mean   :21.89  
          3rd Qu.:30.00  
          Max.   :40.00  

and the dredge() table you used to pick your bestMod1 is in dredgeOut1

dredgeOut1
Global model call: glm(formula = Resp1 ~ Cat1, family = gaussian(link = "identity"), 
    data = myDat1)
---
Model selection table 
  (Intrc) Cat1    R^2 df   logLik  AICc delta weight
2   22.50    + 0.2077  4 -354.584 717.6  0.00      1
1   21.89      0.0000  2 -366.223 736.6 18.98      0
Models ranked by AICc(x) 

11.1.1 Reporting your statistical modelling results - Example 1

Recall that reporting your statistical modelling results means:

  • Reporting your best-specified model(s)

  • Reporting your modelled effects

  • Reporting how well your model explains your response

  • Reporting how important each of your covariates is in explaining your response

Reporting best-specified model(s) - Example 1

Reporting your best-specified model means reporting the terms - the covariates and any interactions - that are in your model.

For Example 1, you will report that your best-specified model is Resp1 ~ Cat1, i.e. that there is evidence that Cat1 explains variability in Resp1. Remember from last week that you can also use the output from dredgeOut1 to report evidence for how you picked the best-specified model, e.g. if you want to report more than one model.

Global model call: glm(formula = Resp1 ~ Cat1, family = gaussian(link = "identity"), 
    data = myDat1)
---
Model selection table 
  (Intrc) Cat1    R^2 df   logLik  AICc delta weight
2   22.50    + 0.2077  4 -354.584 717.6  0.00      1
1   21.89      0.0000  2 -366.223 736.6 18.98      0
Models ranked by AICc(x) 

Reporting your modelled effects - Example 1

Reporting the coefficients

If you remember from last week, you reported your modelled effects by reporting your model coefficients as your coefficients tell you about the effect your covariate has on your response.

With a continuous covariate, the coefficient is the slope of the line (on the linked scale) showing the amount of change in the response that is caused by a change in your continuous covariate.

With a categorical covariate, the coefficient is an intercept, or the model prediction of the response at that level of the categorical covariate (e.g. the modelled Resp1 when Cat1 is Sp1, Sp2 and Sp3). For categorical covariates, there is a coefficient estimated for each level of the covariate that can be thought of as an intercept. So as Cat1 in Example 1 has three levels (Sp1, Sp2 and Sp3), there will be three coefficients estimated (one for each).

Here we will focus on three things:

  • What are your coefficient estimates (with uncertainty)?

  • Is each coefficient estimate different from zero?

  • Are the coefficient estimates different than one another across categorical levels?

What are your coefficient estimates (with uncertainty)?

Last week, you found your coefficients in the “Estimate” column of the summary() output of your model:

coef(summary(bestMod1)) # extract the coefficients from summary()
             Estimate Std. Error   t value     Pr(>|t|)
(Intercept) 22.500000   1.609663 13.978083 5.653157e-25
Cat1Sp2      4.194444   2.146217  1.954343 5.354014e-02
Cat1Sp3     -5.888889   2.146217 -2.743846 7.234664e-03

Interpreting the coefficients from this output takes practice - especially for categorical covariates because of the way that R treats categorical covariates in regression. With R’s “dummy coding”, one level of the covariate (here Sp1) is incorporated into the intercept estimate 22.5 as the reference level. The other coefficients in the Estimate column represent the change in modelled response when you move from the reference level (Sp1) to another level.

The modelled response when Cat1 = Sp1 is 22.5 units.

The modelled response when Cat1 = Sp2 is 4.2 units higher than this reference level = 22.5 + 4.2 = 26.7 units.

The modelled Resp1 when Cat1 = Sp3 is -5.9 units lower than the reference level = 22.5 + -5.9 = 16.6 units.48

So all the information we need is in this summary() output, but not easy to see immediately. An easier way is to use the emmeans package which helps us by reporting the coefficients directly. “emmeans” is the estimated marginal means. The estimated marginal mean is the model prediction for the response when the covariates are held at particular values. If the value of the covariate is not defined, the response prediction is the average response prediction over all values in that covariate.

In our case, we use the emmeans package to get the mean value of the response at each level of the covariate. For categorical covariates, we do this with the emmeans() function:

library(emmeans) # load the emmeans package

emmeans(object = bestMod1, # your model
        specs = ~ Cat1, # your covariates
        type = "response") # report coefficients on the response scale
 Cat1 emmean   SE df lower.CL upper.CL
 Sp1    22.5 1.61 97     19.3     25.7
 Sp2    26.7 1.42 97     23.9     29.5
 Sp3    16.6 1.42 97     13.8     19.4

Confidence level used: 0.95 

So now we have a modelled value of our response for each level of our categorical covariate.
When Cat1 is Sp1, Resp1 is estimated to be 22.5 ± 1.6 units.
When Cat1 is Sp2, Resp1 is estimated to be 26.7 ± 1.4 units.
When Cat1 is Sp3, Resp1 is estimated to be 16.6 ± 1.4 units.

Note that this is the same information in the summary() output just easier to read.49

Note also that two types of uncertainty are measured here. SE stands for the standard error around the prediction, and is a measure of uncertainty of the average modelled effect. The lower.CL and upper.CL represent the 95% confidence limits of the prediction - so if I observed a new Resp1 at a particular Cat1, there would be a 95% chance it would lie between the bounds of the confidence limits.

Is each coefficient estimate different from zero?

Cat1 is in your best-specified model, so you know that there is evidence it has an effect50 on your response (Resp1). But is there an effect on Resp1 for all levels of Cat1 - are any of the coefficients statistically similar to 0?

forComp <- emmeans(object = bestMod1, # your model
                   specs = pairwise ~ Cat1, # your covariates, with the request for a pairwise test
                   type = "response") # report coefficients on the response scale

test(forComp)$emmeans # shows test of coefficients = 0
 Cat1 emmean   SE df t.ratio p.value
 Sp1    22.5 1.61 97  13.978  <.0001
 Sp2    26.7 1.42 97  18.804  <.0001
 Sp3    16.6 1.42 97  11.701  <.0001

Here we see that the coefficients for all levels are significantly different than zero (t-test, P < 0.0001.)

Are the coefficient estimates different than one another across categorical levels?

Finally, we ask if all levels of Cat1 affect Resp1 in the same way - i.e. are the coefficients across factor levels significantly different from one another?

To get evidence about how each level in your categorical covariate affects your response, you need to test which effects differ from one another using multiple comparisons51, i.e. you will compare the modelled effect of each level of your categorical covariate vs. each other level of your categorical covariate to determine which are different from each other (called pairwise testing). You will do this by testing the null hypothesis that the modelled effects are similar to one another, typically rejecting the null hypothesis when p < 0.05. Remember that the p-value is the probability that we observe a value at least as big as the one we observed even though our null hypothesis is true. In this case, we are looking at the difference between coefficients estimated for two levels of our covariate. The p-value is the probability of getting a difference at least as big as the one we observed even though there is actually no difference between the coefficients (the null hypothesis is true).

A couple of things to note about multiple comparison testing:

  1. Multiple comparison testing on a categorical covariate should only be done after your hypothesis testing has given you evidence that the covariate has an effect on your response. That is, you should only do a multiple comparison test on a covariate if that covariate is in your best-specified model. As this is a test done after your hypothesis testing, it is called a post-hoc test.

  2. Multiple comparison testing can be a problem because you are essentially repeating a hypothesis test many times on the same data (i.e. are the effects of Sp1 different than those of Sp2? are the effects of Sp1 different than those of Sp3? are the effects of Sp2 different than those of Sp3?…). These repeated tests mean there is a high chance that you will find a difference in one test purely due to random chance, vs. due to there being an actual difference. To account for this, the multiple comparison tests you will perform have been formulated to correct for this increased error.

Multiple comparison testing is very simple with the emmeans package. It just requires us to add that we would like to have multiple comparisons testing in our call by adding pairwise to the specs = argument:

forComp <- emmeans(object = bestMod1, # your model
                   specs = pairwise ~ Cat1, # your covariates, with the request for a pairwise test
                   type = "response") # report coefficients on the response scale

test(forComp)
$emmeans
 Cat1 emmean   SE df t.ratio p.value
 Sp1    22.5 1.61 97  13.978  <.0001
 Sp2    26.7 1.42 97  18.804  <.0001
 Sp3    16.6 1.42 97  11.701  <.0001


$contrasts
 contrast  estimate   SE df t.ratio p.value
 Sp1 - Sp2    -4.19 2.15 97  -1.954  0.1292
 Sp1 - Sp3     5.89 2.15 97   2.744  0.0196
 Sp2 - Sp3    10.08 2.01 97   5.023  <.0001

P value adjustment: tukey method for comparing a family of 3 estimates 

You can see that you get two tables from this new emmeans() call. The first table ($emmeans) is the same data you saw before - the modelled coefficients (intercepts) for each level in Cat1. The second table ($contrasts) shows the results of the multiple comparison (pairwise) testing. The values in the p.value column tell us the results of the hypothesis test comparing the coefficients between the two levels. For example, for Sp1 vs. Sp2, there is a 13% (p = 0.13) probability of getting a difference in coefficients at least as big as 22.5 - 26.7 = -4.2 even though the null hypothesis (no difference) is true. This value 13% (p = 0.13) is too big (i.e. bigger than our threshold of 5% or p = 0.05) for us to believe we have evidence the coefficients are different from one another.

Based on a threshold p-value of 0.05, we can see that:

The value of Resp1 when Cat1 is Sp1 is not statistically different than that when Cat1 is Sp2 as p = 0.13 is greater than p = 0.05.
The value of Resp1 when Cat1 is Sp1 is different (higher) than that when Cat1 is Sp3 as p = 0.02 is smaller than p = 0.05. The value of Resp1 when Cat1 is Sp2 is different (higher) than that when Cat1 is Sp3 as p < 0.000152 is smaller than p = 0.05.53.

Note that the p-values have been adjusted via the Tukey method which adjusts the difference that the two coefficients need to have to allow for the fact that we’re making multiple comparisons.54

Note that you can also get the results from emmeans visually with

plot(forComp,
     comparisons = TRUE)

Plotting modelled effects

And this brings us to another important aspect of reporting: visualizing your modelled effects. As mentioned last class, you want to include i) your model fit, ii) uncertainty around that fit, and iii) your observations on your plot. You can do this “by hand” with:

# Set up your covariates for the visualized fit
forCat1<-unique(myDat1$Cat1) # every value of your categorical covariate

# create a data frame with your covariates
forVis<-expand.grid(Cat1=forCat1) # expand.grid() function makes sure you have all combinations of covariates

# Get your model fit estimates at each value of your covariates
modFit<-predict(bestMod1, # the model
                newdata = forVis, # the covariate values
                type = "response", # make the predictions on the response scale
                se.fit = TRUE) # include uncertainty estimate

forVis$Fit<-modFit$fit # add your fit to the data frame
forVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frame
forVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame

library(ggplot2)

ggplot() + # start ggplot
  geom_point(data = myDat1, # add observations to your plot
             mapping = aes(x = Cat1, y = Resp1), 
             position=position_jitter(width=0.1)) + # control position of data points so they are easier to see on the plot
  geom_errorbar(data = forVis, # add the uncertainty to your plot
              mapping = aes(x = Cat1, y = Fit, ymin = Lower, ymax = Upper),
              size=1.2) + # control thickness of errorbar line
  geom_point(data = forVis, # add the modelled fit to your plot
              mapping = aes(x = Cat1, y = Fit), 
             shape = 22, # shape of point
             size = 3, # size of point
             fill = "white", # fill colour for plot
             col = 'black') + # outline colour for plot
  ylab("Resp1, (units)")+ # change y-axis label
  xlab("Cat1, (units)")+ # change x-axis label
  theme_bw() # change ggplot theme

Reporting how well your model explains your response - Example 1

As a reminder from last week, you can estimate the deviance explained in your response by using a pseudo \(R^2\) called the Likelihood Ratio \(R^2\) or \(R^2_{LR}\).55. You can calculate \(R^2_{LR}\) by hand, read it from our dredge() output, or estimate it using r.squaredLR() from the MuMIn package:

r.squaredLR(bestMod1)
[1] 0.2076829
attr(,"adj.r.squared")
[1] 0.2078199

Note here that two values of \(R^2_{LR}\) are reported. The adjusted pseudo \(R^2\) given here under attr(,"adj.r.squared") has been scaled so that \(R^2_{LR}\) can reach a maximum of 1, to be equivalent to a regular \(R^2\).

So you estimate that ~ 21% of the deviance in Resp1 is explained by Cat1 via bestMod1.

One nice feature of the \(R^2_{LR}\) is that it is equivalent to the regular \(R^2\) when your model assumes a normal error distribution assumption and linear shape assumption, so you can use \(R^2_{LR}\) for any of the models we’re discussing in class. Let’s check this here by comparing our \(R^2_{LR}\) to that from true \(R^2\) we calculated in Week 9:

1-summary(bestMod1)$deviance/summary(bestMod1)$null.deviance
[1] 0.2076829

\(R^2_{LR}\) and \(R^2\) are identical because you have a normal error distribution assumption and linear shape assumption.

Reporting covariate importance - Example 1

With Example 1, you have only one covariate (Cat1) and so this covariate is responsible for explaining all of the variability in your response (Resp1). You can see that it appears in all models with any weight with your sw() function from the MuMIn package.

dredgeOut1
Global model call: glm(formula = Resp1 ~ Cat1, family = gaussian(link = "identity"), 
    data = myDat1)
---
Model selection table 
  (Intrc) Cat1    R^2 df   logLik  AICc delta weight
2   22.50    + 0.2077  4 -354.584 717.6  0.00      1
1   21.89      0.0000  2 -366.223 736.6 18.98      0
Models ranked by AICc(x) 
sw(dredgeOut1)
                     Cat1
Sum of weights:      1   
N containing models: 1   

11.1.2 Predicting with Example 1

With Example 1, we have only a categorical covariate. In such cases, it doesn’t make sense to use your model to make a prediction as you can’t make a prediction for a level not already in your covariate, e.g. there’s no way to predict Resp1 if Cat1 was Sp5 or Sp89.

So if you have only categorical covariates, you aren’t able to use your model for prediction. Coming up (Example 3) we’ll look at an example with a combination of continuous and categorical covariates and, in these cases, prediction is possible.

11.2 Example 2: Resp2 ~ Cat2 + Cat3 + Cat2:Cat3

For Example #2, assume your best-specified model is that your response variable (Resp2) is explained by two categorical covariates (Cat2 & Cat3) as well as the interaction between the covariates:

Resp2 ~ Cat2 + Cat3 + Cat2:Cat3

Your best-specified model for example #2 is in an object called bestMod2:

bestMod2

Call:  glm(formula = Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1, family = gaussian(link = "identity"), 
    data = myDat2)

Coefficients:
          (Intercept)              Cat2TypeB              Cat2TypeC  
               364.91                  37.37                 -75.94  
            Cat2TypeD             Cat3Treat2            Cat3Control  
              -141.97                 -43.63                  64.12  
 Cat2TypeB:Cat3Treat2   Cat2TypeC:Cat3Treat2   Cat2TypeD:Cat3Treat2  
                42.71                  73.61                  66.02  
Cat2TypeB:Cat3Control  Cat2TypeC:Cat3Control  Cat2TypeD:Cat3Control  
                81.25                 -19.08                 -71.72  

Degrees of Freedom: 99 Total (i.e. Null);  88 Residual
Null Deviance:      1249000 
Residual Deviance: 377000   AIC: 1133

that was fit to data in myDat2:

summary(myDat2)
     Resp2          Cat2         Cat3   
 Min.   :119.9   TypeA:27   Treat1 :29  
 1st Qu.:275.6   TypeB:22   Treat2 :19  
 Median :338.8   TypeC:29   Control:52  
 Mean   :348.3   TypeD:22               
 3rd Qu.:420.3                          
 Max.   :621.5                          

and the dredge() table you used to pick your bestMod2 is in dredgeOut2

dredgeOut2
Global model call: glm(formula = Resp2 ~ Cat2 + Cat3 + Cat2:Cat3, family = gaussian(link = "identity"), 
    data = myDat2)
---
Model selection table 
  (Int) Ct2 Ct3 Ct2:Ct3     R^2 df   logLik   AICc delta weight
8 364.9   +   +       + 0.69800 13 -553.634 1137.5  0.00  0.964
4 354.9   +   +         0.62540  7 -564.419 1144.1  6.55  0.036
2 388.5   +             0.55300  5 -573.243 1157.1 19.63  0.000
1 348.3                 0.00000  2 -613.508 1231.1 93.64  0.000
3 331.1       +         0.03933  4 -611.502 1231.4 93.92  0.000
Models ranked by AICc(x) 

11.2.1 Reporting Example 2

Reporting best-specified model(s) - Example 2

For Example 2, you will report that your best-specified model is Resp2 ~ Cat2 + Cat3 + Cat2:Cat3, i.e. that there is evidence that Cat2 and Cat3 explain variability in Resp2, and that there is an interaction between Cat2 and Cat3 - i.e. the effect of Cat2 on Resp2 depends on the value of Cat3. Remember that you can also use the output from dredgeOut2 to report evidence for how you picked the best-specified model, e.g. if you want to report more than one model.

Global model call: glm(formula = Resp2 ~ Cat2 + Cat3 + Cat2:Cat3, family = gaussian(link = "identity"), 
    data = myDat2)
---
Model selection table 
  (Int) Ct2 Ct3 Ct2:Ct3     R^2 df   logLik   AICc delta weight
8 364.9   +   +       + 0.69800 13 -553.634 1137.5  0.00  0.964
4 354.9   +   +         0.62540  7 -564.419 1144.1  6.55  0.036
2 388.5   +             0.55300  5 -573.243 1157.1 19.63  0.000
1 348.3                 0.00000  2 -613.508 1231.1 93.64  0.000
3 331.1       +         0.03933  4 -611.502 1231.4 93.92  0.000
Models ranked by AICc(x) 

Reporting modelled effects - Example 2

Reporting coefficients
What are your coefficient estimates (with uncertainty)?

Again, for categorical covariates, there is a coefficient estimated for each level of the covariate. If there is one or more interactions among covariates, there will be one coefficient for each combination of levels among covariates. Let’s look at the summary() output of your model to understand better:

coef(summary(bestMod2)) # extract the coefficients from summary()
                        Estimate Std. Error    t value     Pr(>|t|)
(Intercept)            364.90879   24.73850 14.7506420 1.629062e-25
Cat2TypeB               37.37325   32.98467  1.1330491 2.602712e-01
Cat2TypeC              -75.94350   33.87459 -2.2419017 2.748208e-02
Cat2TypeD             -141.96558   38.32472 -3.7042819 3.695537e-04
Cat3Treat2             -43.62974   36.41409 -1.1981554 2.340735e-01
Cat3Control             64.11573   30.29835  2.1161457 3.715672e-02
Cat2TypeB:Cat3Treat2    42.70879   53.60009  0.7968045 4.277091e-01
Cat2TypeC:Cat3Treat2    73.60837   54.15228  1.3592849 1.775295e-01
Cat2TypeD:Cat3Treat2    66.02110   55.13228  1.1975037 2.343260e-01
Cat2TypeB:Cat3Control   81.25030   43.24327  1.8789121 6.356678e-02
Cat2TypeC:Cat3Control  -19.07730   41.29748 -0.4619483 6.452584e-01
Cat2TypeD:Cat3Control  -71.72444   46.17117 -1.5534463 1.239057e-01

Comparing this output to that of Example 1 shows many more estimated coefficients for Example 2. This is because we have one coefficient for each level of each covariate, as well as coefficients for each combination of levels of the covariates.

Again, it takes a little practice to read the coefficients from the summary() output. For Example 2:

The modelled prediction for Resp2 when Cat2 is TypeA and Cat3 is Treat1 is 365 units (the intercept).

The modelled prediction for Resp2 when Cat2 is TypeB and Cat3 is Treat1 is 365 + 37 = 402 units.

The modelled prediction for Resp2 when Cat2 is TypeA and Cat3 is Treat2 is 365 - 44 = 321 units.

The modelled prediction for Resp2 when Cat2 is TypeC and Cat3 is Treat2 is 365 - 76 - 44 + 74 = 319 units.

As above, we can use the emmeans package to more easily see these coefficients:

emmeans(object = bestMod2, # your model
        specs = ~ Cat2 + Cat3 + Cat2:Cat3, # your covariates
        type = "response") # report coefficients on the response scale
 Cat2  Cat3    emmean   SE df lower.CL upper.CL
 TypeA Treat1     365 24.7 88      316      414
 TypeB Treat1     402 21.8 88      359      446
 TypeC Treat1     289 23.1 88      243      335
 TypeD Treat1     223 29.3 88      165      281
 TypeA Treat2     321 26.7 88      268      374
 TypeB Treat2     401 32.7 88      336      466
 TypeC Treat2     319 32.7 88      254      384
 TypeD Treat2     245 29.3 88      187      304
 TypeA Control    429 17.5 88      394      464
 TypeB Control    548 21.8 88      504      591
 TypeC Control    334 15.9 88      302      366
 TypeD Control    215 18.9 88      178      253

Confidence level used: 0.95 

So now we have a modelled value of our response for each level of our categorical covariate. For example:

The modelled prediction for Resp2 when Cat2 is TypeA and Cat3 is Treat1 is 365 +/- 25 units.

The modelled prediction for Resp2 when Cat2 is TypeB and Cat3 is Treat1 is 402 +/- 22 units.

The modelled prediction for Resp2 when Cat2 is TypeA and Cat3 is Treat2 is 321 +/- 27 units.

The modelled prediction for Resp2 when Cat2 is TypeC and Cat3 is Treat2 is 319 +/- 33 units.

Is each coefficient estimate different from zero?

Here we find out if there is an effect on Resp1 for all levels of Cat1 - are any of the coefficients statistically similar to 0?

forComp <- emmeans(object = bestMod2, # your model
                   specs = pairwise ~ Cat2 + Cat3 + Cat2:Cat3, # your covariates, with the request for a pairwise test
                   type = "response") # report coefficients on the response scale

test(forComp)$emmeans # shows test of coefficients = 0
 Cat2  Cat3    emmean   SE df t.ratio p.value
 TypeA Treat1     365 24.7 88  14.751  <.0001
 TypeB Treat1     402 21.8 88  18.439  <.0001
 TypeC Treat1     289 23.1 88  12.487  <.0001
 TypeD Treat1     223 29.3 88   7.617  <.0001
 TypeA Treat2     321 26.7 88  12.024  <.0001
 TypeB Treat2     401 32.7 88  12.264  <.0001
 TypeC Treat2     319 32.7 88   9.746  <.0001
 TypeD Treat2     245 29.3 88   8.381  <.0001
 TypeA Control    429 17.5 88  24.526  <.0001
 TypeB Control    548 21.8 88  25.102  <.0001
 TypeC Control    334 15.9 88  21.040  <.0001
 TypeD Control    215 18.9 88  11.397  <.0001

Here we see that the coefficients for all combinations of levels of our two covariates are significantly different than zero (t-test, P < 0.0001.)

Are the coefficient estimates different than one another across categorical levels?

As with Example 1, we can also find out which combinations of covariate levels are leading to statistically different model predictions in Resp2:

forComp <- emmeans(object = bestMod2, # your model
            specs = pairwise ~ Cat2 + Cat3 + Cat2:Cat3, # your covariates, with pairwise to indicate you want multiple comparisons
            type = "response") # report coefficients on the response scale

test(forComp)
$emmeans
 Cat2  Cat3    emmean   SE df t.ratio p.value
 TypeA Treat1     365 24.7 88  14.751  <.0001
 TypeB Treat1     402 21.8 88  18.439  <.0001
 TypeC Treat1     289 23.1 88  12.487  <.0001
 TypeD Treat1     223 29.3 88   7.617  <.0001
 TypeA Treat2     321 26.7 88  12.024  <.0001
 TypeB Treat2     401 32.7 88  12.264  <.0001
 TypeC Treat2     319 32.7 88   9.746  <.0001
 TypeD Treat2     245 29.3 88   8.381  <.0001
 TypeA Control    429 17.5 88  24.526  <.0001
 TypeB Control    548 21.8 88  25.102  <.0001
 TypeC Control    334 15.9 88  21.040  <.0001
 TypeD Control    215 18.9 88  11.397  <.0001


$contrasts
 contrast                      estimate   SE df t.ratio p.value
 TypeA Treat1 - TypeB Treat1    -37.373 33.0 88  -1.133  0.9923
 TypeA Treat1 - TypeC Treat1     75.944 33.9 88   2.242  0.5248
 TypeA Treat1 - TypeD Treat1    141.966 38.3 88   3.704  0.0180
 TypeA Treat1 - TypeA Treat2     43.630 36.4 88   1.198  0.9878
 TypeA Treat1 - TypeB Treat2    -36.452 41.0 88  -0.889  0.9991
 TypeA Treat1 - TypeC Treat2     45.965 41.0 88   1.120  0.9929
 TypeA Treat1 - TypeD Treat2    119.574 38.3 88   3.120  0.0940
 TypeA Treat1 - TypeA Control   -64.116 30.3 88  -2.116  0.6129
 TypeA Treat1 - TypeB Control  -182.739 33.0 88  -5.540  <.0001
 TypeA Treat1 - TypeC Control    30.905 29.4 88   1.051  0.9959
 TypeA Treat1 - TypeD Control   149.574 31.1 88   4.805  0.0004
 TypeB Treat1 - TypeC Treat1    113.317 31.8 88   3.563  0.0277
 TypeB Treat1 - TypeD Treat1    179.339 36.5 88   4.912  0.0002
 TypeB Treat1 - TypeA Treat2     81.003 34.5 88   2.348  0.4517
 TypeB Treat1 - TypeB Treat2      0.921 39.3 88   0.023  1.0000
 TypeB Treat1 - TypeC Treat2     83.338 39.3 88   2.119  0.6110
 TypeB Treat1 - TypeD Treat2    156.947 36.5 88   4.299  0.0025
 TypeB Treat1 - TypeA Control   -26.742 28.0 88  -0.956  0.9982
 TypeB Treat1 - TypeB Control  -145.366 30.9 88  -4.711  0.0005
 TypeB Treat1 - TypeC Control    68.278 27.0 88   2.531  0.3355
 TypeB Treat1 - TypeD Control   186.947 28.9 88   6.477  <.0001
 TypeC Treat1 - TypeD Treat1     66.022 37.3 88   1.769  0.8297
 TypeC Treat1 - TypeA Treat2    -32.314 35.3 88  -0.914  0.9988
 TypeC Treat1 - TypeB Treat2   -112.396 40.1 88  -2.804  0.1965
 TypeC Treat1 - TypeC Treat2    -29.979 40.1 88  -0.748  0.9998
 TypeC Treat1 - TypeD Treat2     43.631 37.3 88   1.169  0.9900
 TypeC Treat1 - TypeA Control  -140.059 29.0 88  -4.828  0.0003
 TypeC Treat1 - TypeB Control  -258.683 31.8 88  -8.134  <.0001
 TypeC Treat1 - TypeC Control   -45.038 28.1 88  -1.605  0.9027
 TypeC Treat1 - TypeD Control    73.631 29.9 88   2.465  0.3757
 TypeD Treat1 - TypeA Treat2    -98.336 39.6 88  -2.481  0.3654
 TypeD Treat1 - TypeB Treat2   -178.418 43.9 88  -4.064  0.0056
 TypeD Treat1 - TypeC Treat2    -96.001 43.9 88  -2.186  0.5636
 TypeD Treat1 - TypeD Treat2    -22.391 41.4 88  -0.541  1.0000
 TypeD Treat1 - TypeA Control  -206.081 34.1 88  -6.043  <.0001
 TypeD Treat1 - TypeB Control  -324.705 36.5 88  -8.894  <.0001
 TypeD Treat1 - TypeC Control  -111.061 33.3 88  -3.335  0.0532
 TypeD Treat1 - TypeD Control     7.609 34.8 88   0.218  1.0000
 TypeA Treat2 - TypeB Treat2    -80.082 42.2 88  -1.895  0.7587
 TypeA Treat2 - TypeC Treat2      2.335 42.2 88   0.055  1.0000
 TypeA Treat2 - TypeD Treat2     75.945 39.6 88   1.916  0.7460
 TypeA Treat2 - TypeA Control  -107.746 31.9 88  -3.374  0.0478
 TypeA Treat2 - TypeB Control  -226.369 34.5 88  -6.562  <.0001
 TypeA Treat2 - TypeC Control   -12.725 31.1 88  -0.409  1.0000
 TypeA Treat2 - TypeD Control   105.945 32.7 88   3.237  0.0693
 TypeB Treat2 - TypeC Treat2     82.417 46.3 88   1.781  0.8238
 TypeB Treat2 - TypeD Treat2    156.026 43.9 88   3.554  0.0285
 TypeB Treat2 - TypeA Control   -27.663 37.1 88  -0.745  0.9998
 TypeB Treat2 - TypeB Control  -146.287 39.3 88  -3.719  0.0172
 TypeB Treat2 - TypeC Control    67.357 36.4 88   1.852  0.7846
 TypeB Treat2 - TypeD Control   186.027 37.8 88   4.923  0.0002
 TypeC Treat2 - TypeD Treat2     73.609 43.9 88   1.677  0.8739
 TypeC Treat2 - TypeA Control  -110.081 37.1 88  -2.967  0.1366
 TypeC Treat2 - TypeB Control  -228.704 39.3 88  -5.815  <.0001
 TypeC Treat2 - TypeC Control   -15.060 36.4 88  -0.414  1.0000
 TypeC Treat2 - TypeD Control   103.609 37.8 88   2.742  0.2240
 TypeD Treat2 - TypeA Control  -183.690 34.1 88  -5.387  <.0001
 TypeD Treat2 - TypeB Control  -302.313 36.5 88  -8.281  <.0001
 TypeD Treat2 - TypeC Control   -88.669 33.3 88  -2.663  0.2624
 TypeD Treat2 - TypeD Control    30.000 34.8 88   0.861  0.9993
 TypeA Control - TypeB Control -118.624 28.0 88  -4.242  0.0030
 TypeA Control - TypeC Control   95.021 23.6 88   4.023  0.0064
 TypeA Control - TypeD Control  213.690 25.7 88   8.299  <.0001
 TypeB Control - TypeC Control  213.644 27.0 88   7.918  <.0001
 TypeB Control - TypeD Control  332.314 28.9 88  11.514  <.0001
 TypeC Control - TypeD Control  118.669 24.7 88   4.809  0.0004

P value adjustment: tukey method for comparing a family of 12 estimates 

Again, the second table ($contrasts) shows the results of the multiple comparison testing. Based on a threshold p-value of 0.05, we can see that the effects at some combinations of our covariates are not statistically different from each other.

For example, the coefficient (predicted Resp2) when Cat3 = Treat1 and Cat2 = TypeA is 365 and this is 37 lower than the coefficient when Cat3 = Treat1 and Cat2 = TypeB, but this difference is not statistically different as p = 0.99 for this comparison.

On the other hand, some combinations of our covariates are statistically different from each other. For example, A comparison of modelled Resp2 when Cat2 = TypeD and Cat3 = Treat2 (245) vs. Cat2 = TypeB and Cat3 = Control (548) shows that they differ by 302 and are statistically different from one another with p < 0.001.

You can also get the results from emmeans visually with

plot(forComp,
     comparisons = TRUE)

Plotting modelled effects

You can visualize your results (model fit, uncertainty and observations) for Example 2 with:

# Set up your covariates for the visualized fit
forCat2<-unique(myDat2$Cat2) # every level of your categorical covariate
forCat3<-unique(myDat2$Cat3) # every level of your categorical covariate
  
# create a data frame with your covariates
forVis<-expand.grid(Cat2 = forCat2, Cat3 = forCat3) # expand.grid() function makes sure you have all combinations of covariates

# Get your model fit estimates at each value of your covariates
modFit<-predict(bestMod2, # the model
                newdata = forVis, # the covariate values
                type = "response", # make the predictions on the response scale
                se.fit = TRUE) # include uncertainty estimate

forVis$Fit<-modFit$fit # add your fit to the data frame
forVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frame
forVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame

library(ggplot2) # load ggplot2 library

ggplot() + # start ggplot
  geom_point(data = myDat2, # add observations to your plot
             mapping = aes(x = Cat2, y = Resp2, col = Cat3), 
             position=position_jitterdodge(jitter.width=0.75, dodge.width=0.75)) + # control position of data points so they are easier to see on the plot
  geom_errorbar(data = forVis, # add the uncertainty to your plot
              mapping = aes(x = Cat2, y = Fit, ymin = Lower, ymax = Upper, col = Cat3),
              position=position_dodge(width=0.75), # control position of data points so they are easier to see on the plot
              size=1.2) + # control thickness of errorbar line
  geom_point(data = forVis, # add the modelled fit to your plot
              mapping = aes(x = Cat2, y = Fit, fill = Cat3), 
             shape = 22, # shape of point
             size=3, # size of point
             col = 'black', # outline colour for point
             position=position_dodge(width=0.75)) + # control position of data points so they are easier to see on the plot
  ylab("Resp2, (units)")+ # change y-axis label
  xlab("Cat2, (units)")+ # change x-axis label
  labs(fill="Cat3, (units)", col="Cat3, (units)") + # change legend title
  theme_bw() # change ggplot theme

Reporting how well your model explains your response - Example 2

As for Example 1, you can report the explained deviance as the Likelihood Ratio \(R^2\):

r.squaredLR(bestMod2)
[1] 0.6980474
attr(,"adj.r.squared")
[1] 0.6980507

So ~ 70% of the deviance in Resp2 is explained by Cat2 and Cat3 and the interaction between the two covariates.

Reporting covariate importance - Example 2

You can report the importance of each model term in explaining deviance in Resp2 with the sum of Akaike weights function (sw() in the MuMIn package), as we did in class:

dredgeOut2
Global model call: glm(formula = Resp2 ~ Cat2 + Cat3 + Cat2:Cat3, family = gaussian(link = "identity"), 
    data = myDat2)
---
Model selection table 
  (Int) Ct2 Ct3 Ct2:Ct3     R^2 df   logLik   AICc delta weight
8 364.9   +   +       + 0.69800 13 -553.634 1137.5  0.00  0.964
4 354.9   +   +         0.62540  7 -564.419 1144.1  6.55  0.036
2 388.5   +             0.55300  5 -573.243 1157.1 19.63  0.000
1 348.3                 0.00000  2 -613.508 1231.1 93.64  0.000
3 331.1       +         0.03933  4 -611.502 1231.4 93.92  0.000
Models ranked by AICc(x) 
sw(dredgeOut2)
                     Cat2 Cat3 Cat2:Cat3
Sum of weights:      1.00 1.00 0.96     
N containing models:    3    3    1     

This tells you that Cat2 and Cat3 appear in all models with any weight and that the interaction term Cat2:Cat3 is slightly less important in explaining deviance in Resp2, though it does appear in the best-specified model.

11.2.2 Predicting with Example 2

As with Example 1, it doesn’t make sense to make predictions with your model for Example 2 as you can’t make predictions to levels of your categorical covariates not already in your model.

11.3 Example 3: Resp3 ~ Cat4 + Cont5 + Cat4:Cont5

For Example #3, assume your best-specified model is that your response variable (Resp3) is explained by one categorical covariate (Cat4) and one continuous covariate (Cont5) as well as the interaction between the covariates:

Resp3 ~ Cat4 + Cont5 + Cat4:Cont5

Your best-specified model for example #3 is in an object called bestMod3:

bestMod3

Call:  glm(formula = Resp3 ~ Cat4 + Cont5 + Cat4:Cont5 + 1, family = gaussian(link = "identity"), 
    data = myDat3)

Coefficients:
    (Intercept)        Cat4Urban         Cat4Wild            Cont5  
      98.346793        -0.236424        -0.867336        -0.002931  
Cat4Urban:Cont5   Cat4Wild:Cont5  
       0.015117         0.030858  

Degrees of Freedom: 99 Total (i.e. Null);  94 Residual
Null Deviance:      4454 
Residual Deviance: 642.8    AIC: 483.9

that was fit to data in myDat3:

summary(myDat3)
     Resp3           Cat4        Cont5      
 Min.   : 92.19   Farm :34   Min.   :310.9  
 1st Qu.: 97.88   Urban:39   1st Qu.:418.7  
 Median :104.22   Wild :27   Median :497.8  
 Mean   :103.98              Mean   :509.8  
 3rd Qu.:108.63              3rd Qu.:609.6  
 Max.   :118.41              Max.   :695.5  

and the dredge() table you used to pick your bestMod3 is in dredgeOut3

dredgeOut3
Global model call: glm(formula = Resp3 ~ Cat4 + Cont5 + Cat4:Cont5, family = gaussian(link = "identity"), 
    data = myDat3)
---
Model selection table 
   (Int) Ct4       Cn5 Ct4:Cn5    R^2 df   logLik  AICc  delta weight
8  98.35   + -0.002931       + 0.8557  7 -234.930 485.1   0.00      1
4  92.25   +  0.009650         0.8166  5 -246.915 504.5  19.39      0
2  96.93   +                   0.7923  4 -253.133 514.7  29.61      0
3  94.92      0.017780         0.0848  3 -327.278 660.8 175.73      0
1 104.00                       0.0000  2 -331.708 667.5 182.46      0
Models ranked by AICc(x) 

11.3.1 Reporting Example 3

Reporting best-specified model(s) - Example 3

For Example 3, you will report that your best-specified model is Resp3 ~ Cat4 + Cont5 + Cat4:Cont5. This model says that there is evidence that Cat4 and Cont5 explain variability in Resp3, and that there is an interaction between Cat4 and Cont5 - i.e. the effect of Cont5 on Resp3 depends on the value of Cat4. Remember that you can also use the output from dredgeOut3 to report evidence for how you picked the best-specified model, e.g. if you want to report more than one model.

Global model call: glm(formula = Resp3 ~ Cat4 + Cont5 + Cat4:Cont5, family = gaussian(link = "identity"), 
    data = myDat3)
---
Model selection table 
   (Int) Ct4       Cn5 Ct4:Cn5    R^2 df   logLik  AICc  delta weight
8  98.35   + -0.002931       + 0.8557  7 -234.930 485.1   0.00      1
4  92.25   +  0.009650         0.8166  5 -246.915 504.5  19.39      0
2  96.93   +                   0.7923  4 -253.133 514.7  29.61      0
3  94.92      0.017780         0.0848  3 -327.278 660.8 175.73      0
1 104.00                       0.0000  2 -331.708 667.5 182.46      0
Models ranked by AICc(x) 

Reporting modelled effects - Example 3

Reporting coefficients

The coefficients estimated for this model will include values of the response (Resp3) at each level of the categorical covariate (Cat4) as well as slopes describing the change in the response when the continuous covariate (Cont5) changes, and a different slope will be estimated for each level in Cat4 (this is because of the interaction in the model).

What are your coefficient estimates (with uncertainty)?

As above, you can take a look at your modelled coefficients with the summary() output:

coef(summary(bestMod3))
                   Estimate  Std. Error     t value     Pr(>|t|)
(Intercept)     98.34679337 1.868900753 52.62280151 1.536692e-71
Cat4Urban       -0.23642403 2.895404981 -0.08165491 9.350947e-01
Cat4Wild        -0.86733554 3.238562391 -0.26781498 7.894285e-01
Cont5           -0.00293078 0.003740830 -0.78345704 4.353283e-01
Cat4Urban:Cont5  0.01511745 0.005611511  2.69400763 8.361918e-03
Cat4Wild:Cont5   0.03085821 0.006183238  4.99062262 2.756559e-06

This shows:

The model prediction of Resp3 when Cat4 is Farm and Cont5 is 0 is 98.34 units56 (the intercept).

The model prediction of Resp3 when Cat4 is Urban and Cont5 is 0 is 98.34 - 0.24 = 98.1 units.

The model prediction of Resp3 when Cat4 is Wild and Cont5 is 0 is 98.34 - 0.87 = 97.5 units.

The slope of the relationship between Cont5 and Resp3 when Cat4 is Farm is -0.0029.57

The slope of the relationship between Cont5 and Resp3 when Cat4 is Urban is -0.0029 + 0.015 = 0.0121.

The slope of the relationship between Cont5 and Resp3 when Cat4 is Wild is -0.0029 + 0.031 = 0.0281.

Again, interpreting the coefficients from the summary() output is tedious and not necessary: You can use the emmeans package to give you the modelled response for each level of the categorical covariate (Cat4) directly:

emmeans(object = bestMod3, # your model
        specs = ~ Cat4 + Cont5 + Cat4:Cont5, # your covariates
        type = "response") # report coefficients on the response scale
 Cat4  Cont5 emmean    SE df lower.CL upper.CL
 Farm    510   96.9 0.458 94     95.9     97.8
 Urban   510  104.3 0.421 94    103.5    105.2
 Wild    510  111.7 0.511 94    110.7    112.7

Confidence level used: 0.95 

Note that emmeans() sets our continuous covariate (Cont5) to the mean value of Cont5 (510 units). We can also control this with the at = argument:

emmeans(object = bestMod3, # your model
        specs = ~ Cat4 + Cont5 + Cat4:Cont5, # your covariates
        type = "response",  # report coefficients on the response scale
        at = list(Cont5 = 0)) # control the value of your continuous covariate at which to make the coefficient estimates
 Cat4  Cont5 emmean   SE df lower.CL upper.CL
 Farm      0   98.3 1.87 94     94.6      102
 Urban     0   98.1 2.21 94     93.7      103
 Wild      0   97.5 2.64 94     92.2      103

Confidence level used: 0.95 

By setting at = 0, you get the intercept - i.e. the modelled Resp3 when Cont5 = 0 for each level of Cat4, and this is what is reported in the summary() output.

Similarly, you can get the emmeans package to give you the slope coefficients for the continuous covariate (Cont5) using the emtrends() function, rather than interpretting them from the summary() output:

emtrends(bestMod,  ~ Cat4, # your categorical covariate
         var = "Cont5") # your continuous covariates
 Cat4  Cont5.trend      SE df lower.CL upper.CL
 Farm     -0.00293 0.00374 94 -0.01036   0.0045
 Urban     0.01219 0.00418 94  0.00388   0.0205
 Wild      0.02793 0.00492 94  0.01815   0.0377

Confidence level used: 0.95 

One thing to note:

The emmeans() function will convert the coefficients (intercepts) from the link to the response scale. You can ask for this (as we did above) with the type = "response" argument. Note that it makes no difference for these examples as we are using a normal error distribution assumption and so the link and response scale are identical (i.e. when link = "identity").

In contrast, the emtrends() function does not convert the coefficients (slopes) to represent effects on the response scale. This is because emtrends() is reporting the slope of a straight line - the trend line on the link scale. But the line isn’t straight on the response scale.


We need to convert from the link to the response by hand when we use emtrends()^[remember, the coefficients for categorical covariates are also being converted from the link to the response scale. It is just that R is doing it for us in the emmeans() function when we add the argument `type = “response”). How we make the conversion, and interpret the result, will depend on our error distribution assumption:

Again, if we are using link = "identity" (as with a normal error distribution), the link and response scale are identical and no conversion is necessary. In this case, the coefficient tells you the absolute change in your response from a unit change in your covariate.

If you are using link = "log" (as with a poisson error distribution, and sometimes also a Gamma error distribution), you get your coefficient on the response scale by taking e to the coefficient on the link scale (with the R function exp()). This coefficient is called the rate ratio58 and it tells you the % change in the response for a unit change in your covariate.

But why do we convert coefficients from the link to response scale with ex when link = "log"?

Given

\[ \begin{align} log_e(\mu_i|Cov_i) &= \beta_1 \cdot Cov_i + \beta_0 \\ Resp_i &\sim poisson(\mu_i) \end{align} \]

Then the rate ratio (RR), or % change in the response for a unit change in the covariate becomes,

\[ \begin{align} RR&=\frac{(\mu_i|Cov = a+1)}{(\mu_i|Cov = a)}\\[2em] log_e(RR)&=log_e(\frac{(\mu_i|Cov = a+1)}{(\mu_i|Cov = a)})\\[2em] log_e(RR)&=log_e(\mu_i|Cov = a+1)-log_e(\mu_i|Cov = a)\\[2em] log_e(RR)&=(\beta_1\cdot (a+1)+\beta_0)-(\beta_1\cdot (a)+\beta_0)\\[2em] log_e(RR)&=\beta_1\cdot a + \beta_1+\beta_0-\beta_1\cdot a+\beta_0\\[2em] log_e(RR)&=\beta_1\\[2em] RR &=exp^{\beta_1} \end{align} \]

If you are using link = "logit" (as with a binomial error distribution), you also get your coefficient on the response scale by taking e to the coefficient on the link scale (with the R function exp()), but here your coefficient on the response scale tells you the % change in the odds (probability of a success vs. the probability of a failure) given a unit change in your covariate - this estimate is called the odds ratio.

But why do we convert coefficients from the link to response scale with ex when link = "logit"?

Given a binomial model, \[ \begin{align} log_e(\frac{p_i}{1-p_i}) &= \beta_1 \cdot Cov_i + \beta_0 \\ Resp_i &\sim binomial(p_i) \end{align} \] \(p_i\) is the probability of success and \(1-p_i\) is the probability of failure, and the odds are \(\frac{p_i}{1-p_i}\).

Then the odds ratio (OR), or % change in odds due to a change in covariate is:

\[ \begin{align} OR&=\frac{odds\;when\;Cov=a+1}{odds\;when \;Cov=a}\\[2em] OR&=\frac{(\frac{p}{1-p}|Cov=a+1)}{(\frac{p}{1-p}|Cov=a)}\\[2em] log_e(OR)&=(log_e\frac{p}{1-p}|Cov=a+1) - (log_e\frac{p}{1-p}|Cov=a)\\[2em] log_e(OR)&=(\beta_1\cdot (a+1)+\beta_0) - (\beta_1\cdot a+\beta_0)\\[2em] log_e(OR)&=\beta_1\cdot a + \beta_1-\beta_1\cdot a\\[2em] OR &=exp^{\beta_1} \end{align} \]


For a model where link = "inverse", the interpretation on the response scale is less clear. In this case, you would use your model to make predictions and report the effects by describing these changes in predictions.

Is each coefficient estimate different from zero?

You can find out if the coefficients are different than zero, similar to above. You use emmeans() with test() when the covariate is categorical;

forInt <- emmeans(object = bestMod3, # your model
                  specs = ~ Cat4 + Cont5 + Cat4:Cont5, # your covariates
                  type = "response",  # report coefficients on the response scale
                  at = list(Cont5 = 0)) # control the value of your continuous covariate at which to make the coefficient estimates

test(forInt) # get test if coefficient is different than zero.
 Cat4  Cont5 emmean   SE df t.ratio p.value
 Farm      0   98.3 1.87 94  52.623  <.0001
 Urban     0   98.1 2.21 94  44.364  <.0001
 Wild      0   97.5 2.64 94  36.856  <.0001

and we see that all coefficients are significantly different than zero (t-test, P < 0.0001).

For the continuous covariate, you use emtrends() with test():

forSlope <- emtrends(bestMod, # your model
                      ~ Cat4, # your categorical covariate, with a request for pairwise testing
                      var = "Cont5") # your continuous covariate

test(forSlope)
 Cat4  Cont5.trend      SE df t.ratio p.value
 Farm     -0.00293 0.00374 94  -0.783  0.4353
 Urban     0.01219 0.00418 94   2.914  0.0045
 Wild      0.02793 0.00492 94   5.673  <.0001

And we can see that

  • the slope associated with Cont5 when Cat4 = Farm is not different than 0 (p = 0.44). This means that there is no effect of Cont5 on Resp3 when Cat4 = Farm.

  • the slope associated with Cont5 when Cat4 = Urban is different than 0 (p = 0.0045) and positive (0.012). This means that Resp3 increases with Cont5 when Cat4 = Urban.

  • the slope associated with Cont5 when Cat4 = Wild is different than 0 (p < 0.0001) and positive (0.028). This means that Resp3 increases with Cont5 when Cat4 = Wild.

Are the coefficient estimates different than one another across categorical levels?

And as with the other examples, you can find which coefficients are significantly different from eachother across the levels of your categorical covariate. A rule here is to always check for differences in the coefficients associated with your continuous covariate(s) first as these differences are slope differences. If you have a difference in slope it is very likely you have intercept differences because the lines are not parallel.


You can also find out which slopes (i.e. the effect of Cont5 on Resp3) are different across the levels of Cat4 using emtrends() with a request for pairwise testing:

forCompSlope <- emtrends(bestMod, # your model
                         pairwise ~ Cat4, # your categorical covariate, with a request for pairwise testing
                         var = "Cont5") # your continuous covariate

forCompSlope
$emtrends
 Cat4  Cont5.trend      SE df lower.CL upper.CL
 Farm     -0.00293 0.00374 94 -0.01036   0.0045
 Urban     0.01219 0.00418 94  0.00388   0.0205
 Wild      0.02793 0.00492 94  0.01815   0.0377

Confidence level used: 0.95 

$contrasts
 contrast     estimate      SE df t.ratio p.value
 Farm - Urban  -0.0151 0.00561 94  -2.694  0.0226
 Farm - Wild   -0.0309 0.00618 94  -4.991  <.0001
 Urban - Wild  -0.0157 0.00646 94  -2.437  0.0437

P value adjustment: tukey method for comparing a family of 3 estimates 

From $contrasts you see that all slope estimates for all levels of Cat4 are significantly different from one another (0.0001 ≤ p ≤ 0.044).

Note that you can visualize these differences in slopes with:

plot(forCompSlope,
     comparisons = TRUE)

If slopes were similar, you will want to test if the levels of your categorical covariate (Cat4) have significantly different coefficients (intercepts) from each other with emmeans():

forCompInt <- emmeans(object = bestMod3, # your model
            specs = pairwise ~ Cat4 + Cont5 + Cat4:Cont5, # your covariates, with the request for a pairwise test
            type = "response") # report coefficients on the response scale

forCompInt
$emmeans
 Cat4  Cont5 emmean    SE df lower.CL upper.CL
 Farm    510   96.9 0.458 94     95.9     97.8
 Urban   510  104.3 0.421 94    103.5    105.2
 Wild    510  111.7 0.511 94    110.7    112.7

Confidence level used: 0.95 

$contrasts
 contrast                               estimate    SE df t.ratio p.value
 Farm Cont5509.767 - Urban Cont5509.767    -7.47 0.622 94 -12.014  <.0001
 Farm Cont5509.767 - Wild Cont5509.767    -14.86 0.686 94 -21.667  <.0001
 Urban Cont5509.767 - Wild Cont5509.767    -7.39 0.662 94 -11.175  <.0001

P value adjustment: tukey method for comparing a family of 3 estimates 
plot(forCompInt,
     comparisons = TRUE)

From the bottom table ($contrasts), you can see that the intercept (modelled Resp2) modelled for each level in Cat4 is significantly different than every other level (all p-values are < 0.0001).

Plotting modelled effects

You can visualize your results (model fit, uncertainty and observations) for Example 3 with (the code is different than the other examples as you have a continuous covariate now):

# Set up your covariates for the visualized fit
forCat4<-unique(myDat3$Cat4) # every level of your categorical covariate
forCont5<-seq(from = min(myDat3$Cont5), to = max(myDat3$Cont5), by = 1)# a sequence of numbers in your continuous covariate range
  
# create a data frame with your covariates
forVis<-expand.grid(Cat4 = forCat4, Cont5 = forCont5) # expand.grid() function makes sure you have all combinations of covariates

# Get your model fit estimates at each value of your covariates
modFit<-predict(bestMod3, # the model
                newdata = forVis, # the covariate values
                type = "response", # make the predictions on the response scale
                se.fit = TRUE) # include uncertainty estimate

forVis$Fit<-modFit$fit # add your fit to the data frame
forVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frame
forVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame

library(ggplot2) # load ggplot2 library

ggplot() + # start ggplot
  
  geom_point(data = myDat3, # add observations to your plot
             mapping = aes(x = Cont5, y = Resp3, col = Cat4)) + # control position of data points so they are easier to see on the plot
  
  geom_line(data = forVis, # add the modelled fit to your plot
              mapping = aes(x = Cont5, y = Fit, col = Cat4),
              size = 1.2) + # control thickness of line
  
    geom_line(data = forVis, # add uncertainty to your plot (upper line)
              mapping = aes(x = Cont5, y = Upper, col = Cat4),
              size = 0.8, # control thickness of line
              linetype = 2) + # control style of line
  
      geom_line(data = forVis, # add uncertainty to your plot (lower line)
              mapping = aes(x = Cont5, y = Lower, col = Cat4),
              size = 0.8, # control thickness of line
              linetype = 2) + # control style of line
  
  ylab("Resp3, (units)") + # change y-axis label
  
  xlab("Cont5, (units)") + # change x-axis label
  
  labs(col="Cat4, (units)") + # change legend title
  
  theme_bw() # change ggplot theme

As you saw in the previous section, the slope describing the effect of Cont5 on Resp3 when Cat4 is Farm was not significantly different than 0 (p = 0.44). Usually you would then not plot this line when plotting effects to avoid misinterpretation, but I’ll leave it on the plot for now.

Reporting how well your model explains your response - Example 3

As with the previous examples, you can report the explained deviance as the Likelihood Ratio \(R^2\):

r.squaredLR(bestMod3)
[1] 0.855657
attr(,"adj.r.squared")
[1] 0.8567834

So 86 % of the deviance in Resp3 is explained by variation in Cat4, Cont5 and the interaction between the two covariates.

Reporting covariate importance - Example 3

As with the previous examples, you can report the importance of each model term in explaining deviance in the response with the sum of Akaike weights function (sw() in the MuMIn package):

dredgeOut3
Global model call: glm(formula = Resp3 ~ Cat4 + Cont5 + Cat4:Cont5, family = gaussian(link = "identity"), 
    data = myDat3)
---
Model selection table 
   (Int) Ct4       Cn5 Ct4:Cn5    R^2 df   logLik  AICc  delta weight
8  98.35   + -0.002931       + 0.8557  7 -234.930 485.1   0.00      1
4  92.25   +  0.009650         0.8166  5 -246.915 504.5  19.39      0
2  96.93   +                   0.7923  4 -253.133 514.7  29.61      0
3  94.92      0.017780         0.0848  3 -327.278 660.8 175.73      0
1 104.00                       0.0000  2 -331.708 667.5 182.46      0
Models ranked by AICc(x) 
sw(dredgeOut3)
                     Cat4 Cont5 Cat4:Cont5
Sum of weights:      1    1     1         
N containing models: 3    3     1         

This tells you that Cat4, Cont5 appear in all models with any weight, indicating all terms have equal importance in explaining the deviance of the response.

11.3.2 Predicting with Example 3

Unlike with the other examples in this week’s notes, with Example 3, the presence of a continuous covariate (Cont5) means you can make predictions with your model - i.e. you can estimate what Resp3 might be if Cont5 was higher or lower. Because you also have a categorical covariate (Cat4) in your model, you do need to specify a level for Cat4 when you make your prediction.

Recall from Week 11 that we used the predict() function to make our predictions. For example, you can estimate the predicted Resp3 when Cont5 = 350 units and Cat4 = Urban with:

predict(object = bestMod, # your model
        newdata = data.frame(Cont5 = 350, Cat4 = "Urban"), # the values of the covariates at which to make the prediction
        type = "response", # to make the prediction on the response scale.
        se.fit = TRUE) # to include a measure of uncertainty around the prediction
$fit
       1 
102.3757 

$se.fit
[1] 0.822138

$residual.scale
[1] 2.6151

So the predicted Resp3 when Cont5 = 350 units and Cat4 = Urban is 102.38 ± 0.82.

To get the predicted Resp3 when Cont5 = 1200 units when Cat4 is Wild, use:

predict(object = bestMod, # your model
        newdata = data.frame(Cont5 = 1200, Cat4 = "Wild"), # the values of the covariates at which to make the prediction
        type = "response", # to make the prediction on the response scale.
        se.fit = TRUE) # to include a measure of uncertainty around the prediction
$fit
       1 
130.9924 

$se.fit
[1] 3.349379

$residual.scale
[1] 2.6151

which shows the predicted Resp3 when Cont5 = 1200 units and Cat4 = Wild is 130.99 ± 3.35.

Footnotes

  1. this is called the Nagelkerke’s modified statistic - see ?r.squaredLR for more information↩︎

  2. see also the course notes from Model Validation↩︎

  3. i.e. in the dredge() output↩︎

  4. you will be using the predict() function again in an upcoming chapter… when you predict!↩︎

  5. note that the p-values reported in the coefficient table are the result of a test of whether the coefficient associated with Sp2 is different than that of Sp1 (p = 0.001), and if the effect of Sp3 is different than that of Sp1 (p = 0.017)), but a comparison of effects of Sp2 vs. Sp3 is missing. We will discuss this more further along in these notes.↩︎

  6. emmeans stands for “estimated marginal means”↩︎

  7. note that the summary() output becomes easier to read if you force the starting model not to have an intercept with: startMod<-glm(formula = Resp1 ~ Cat1 - 1, data = myDat, family = gaussian(link="identity")). This is fine to do here where you only have categorical predictors but causes problems when you start having continuous predictors in your model as well. So we won’t be practicing this in class and instead will use the very flexible and helpful emmeans package to help us report coefficients from our statistical models.↩︎

  8. “post hoc” is a Latin phrase meaning “after the event”↩︎

  9. when the P-value is very low, R reports is as simple <.0001↩︎

  10. note that now you get to assess the difference between coefficients for Sp2 and Sp3 which was missing in the summary() output above↩︎

  11. The emmeans package is very flexible and has a lot of options as to how to make these corrections depending on your needs. Plenty more information is available here↩︎

  12. * in units of Resp3 per units of Cont4↩︎

  13. units of Resp4↩︎

  14. units will be the units of Resp4 divided by those of Cont6↩︎

  15. in units of Resp4 per units of Cont6↩︎

  16. in units of Resp4 per units of Cont6↩︎

  17. in units of Resp4 per units of Cont6↩︎

  18. in units of Resp4 per units of Cont6↩︎

  19. also called the incidence rate ratio, incidence density ratio or relative risk↩︎

  20. a reminder that the number e that the function exp() uses as its base is Euler’s number. It (along with the natural logarithm, log()) is used to model exponential growth or decay, and can be used here when you want models where the response can not be below zero.↩︎

  21. this is what R uses when you say link = "log"↩︎

  22. from the emmeans package↩︎

  23. remember, the coefficients for categorical predictors are also being converted from the link to the response scale. It is just that R is doing it for us in the emmeans() function when we add the argument `type = “response”)↩︎

  24. from the emmeans package↩︎

  25. taken from the emmOut output above↩︎

  26. the mean of the range of Cont6↩︎

  27. the mean of the range of Cont6↩︎

  28. the mean of the range of Cont6↩︎

  29. i.e. effects associated with the continuous predictor↩︎

  30. Notice I write “response(s)” in the title of this section - plural. It is possible to have multiple response variables and we’ll discuss this in class. For the focus of this course though, we’ll be working with one response variable↩︎

  31. here, I’ll write “citation” where you would support your statements with the existing literature↩︎

  32. find this citation with citation("ggplot2")↩︎

  33. you can find which version of R you are using with citation() (nothing inside the parentheses)↩︎

  34. you can find the citation for the car package with citation("car"). Remember to cite every package you use in R.↩︎

  35. see readme file above↩︎

  36. find this with citation("DHARMa")↩︎

  37. As a reminder, in Week 9 we also discussed testing your hypothesis via p-values, as well as the limitations of this method when you have more than one covariate↩︎

  38. AICc is the corrected Akaike Information Criterion which is more conservative than traditional AIC, i.e. models with more covariates need to increase the explained deviance quite a bit before the AICc metric improves↩︎

  39. note that you can also include tables directly from R using the gt package↩︎

  40. watch your number of significant units here. Make sure they make sense for the type of measurement. And be consistent↩︎

  41. here equivalent to the traditional R2 as we have a normal error distribution assumption and linear shape assumption. The traditional R2 is found with 1-summary(bestMod)$deviance/summary(bestMod)$null.deviance↩︎

  42. we’ll go over this in class today, but see also your course notes from week 9↩︎

  43. more to come!↩︎

  44. we will be using the predict() function again… when we predict!↩︎

  45. see [https://plotly.com/r/getting-started/] if you are interested in learning more about the plotly package↩︎

  46. see also the course notes from Week 9↩︎

  47. i.e. in the dredge() output↩︎

  48. note that the p-values reported in the coefficient table are the result of a test of whether the coefficient associated with Sp2 is different than that of Sp1 (p = 0.054), and if the effect of Sp3 is different than that of Sp1 (p = 0.0072), but a comparison of effects of Sp2 vs. Sp3 is missing. We’ll discuss this more further along in these notes.↩︎

  49. note that the summary() output becomes easier to read if you force the starting model not to have an intercept with: startMod<-glm(formula = Resp1 ~ Cat1-1, data = myDat, family = gaussian(link="identity")). This is fine to do here where you only have categorical covariates but causes problems when you start having continuous covariates in your model as well. So we won’t be practicing this in class and instead will use the very flexible and helpful emmeans package to help us report coefficients from our statistical models.↩︎

  50. here I say that your covariate “has an effect” on your response. Your covariate being in your best-specified model just means that there is evidence that variation in your covariate explains variation in your response. It’s up to you to think about the mechanisms underlying this relationship - i.e. is it causation or just correlation?↩︎

  51. see section 11.6 in Crawley for more on multiple comparisons↩︎

  52. when the p-value is very low, R reports is as simple <.0001↩︎

  53. note that now we get to assess the difference between coefficients for Sp2 and Sp3 which was missing in the summary() output above↩︎

  54. The emmeans package is very flexible and has a lot of options as to how to make these corrections depending on your needs. Plenty more information is available here↩︎

  55. see Week 11’s course notes for more on this↩︎

  56. units of Resp3↩︎

  57. units will be the units of Resp3 divided by those of Cont5↩︎

  58. also called the incidence rate ratio, incidence density ratio or relative risk↩︎